Janis Janusevskis
email: janisj@latnet.lv

Rodolphe Le Riche
email: leriche@emse.fr

David Ginsbourger
email: david.ginbourger@stat.unibe.ch

Parallel expected improvements for global optimization: summary, bounds and speed-up

Keywords: Kriging based optimization, Gaussian process, Expected improvement, Optimization using computer grids, Distributed calculations, Monte Carlo optimization

The sequential sampling strategies based on Gaussian processes are widely used for optimization of time consuming simulators. In practice, such computationally demanding problems are solved by increasing number of processing units. This has therefore induced extensions of sampling criteria which consider the framework of parallel calculation.

This report further studies expected improvement criteria for parallel and asynchronous computations. A unified parallel asynchronous expected improvement criterion is formulated. Bounds and strategies for comparing criteria values at various design points are discussed. Finally, the impact of the number of available computing units on the performance is empirically investigated.

Introduction

Many design problems are nowadays based on computationally costly simulator codes. The current technology for adressing computationally intensive tasks is based on an increasing number of processing units (processors, cores, CPUs, GPUs). For example, the next generation of exascale computers (10 18 floating point operations per second) should have of the order of 10 9 computing nodes [START_REF] William | Scientific directions for exascale computing research[END_REF]. Distributing calculations on several nodes is particularly relevant in optimization where a natural parallelization level exists, that of calculating concurrent designs in parallel. However, most current optimization methods have been thought sequentially. Today, there is a need for optimization strategies that operate within the framework of heterogeneous parallel computation, i.e. algorithms designed for being used on computing grids.

The Efficient Global Optimization (EGO) algorithm [START_REF] Donald | A taxonomy of global optimization methods based on response surfaces[END_REF] has become a popular choice for optimizing computationally expensive simulation models. EGO is based on kriging (Gaussian process regression) metamodels, that allow to formulate the Expected improvement (EI) criterion, which provides a compromise between global exploration of the domain space and local exploitation of the best known solutions. At each optimization step the EI is maximized to obtain the next simulation point. By definition the EI is a sequential one point strategy [START_REF] Ginsbourger | Towards GP-based optimization with finite time horizon[END_REF] and lacks means for efficiently accounting for parallel processing capabilities.

The EI criterion has been adapted for parallel computation in [START_REF] Ginsbourger | Kriging is well-suited to parallelize optimization[END_REF], [START_REF] Schonlau | Computer experiments and global optimization[END_REF] leading to the so called "q steps EI" (qEI), that provides the synchronous selection of q new points for the next simulation run. This qEI criterion assumes that the computation time is constant for all design points, i.e. once a set of points is sent for simulation the results are available at the same time for all the points in the set. In general this assumption may not be correct. The simulations at different points may have different algorithmic complexities or run on nodes with different performances or loads, therefore the simulations will finish successively, i.e. new computing nodes become available asynchronously.

The selection of points for the next calculation should be made as soon as computational resources are available, i.e. asynchronously. In [START_REF] Ginsbourger | Dealing with asynchronicity in parallel Gaussian Process based global optimization[END_REF] the qEI criterion has been extended and the so called "Expected Expected Improvement" (EEI) criterion has been proposed. In the framework of Gaussian processes one assumes that the simulator response at busy points (points that have been sent to simulator but have not returned yet) are conditional Gaussian random variables. Maximizing expectation of the EI with respect to busy points allows selection of a new set of points accounting for the points that are being simulated.

As already noted in [START_REF] Ginsbourger | Kriging is well-suited to parallelize optimization[END_REF] and [START_REF] Ginsbourger | Dealing with asynchronicity in parallel Gaussian Process based global optimization[END_REF], the implementation of these criteria faces several problems. Firstly, in general qEI and EEI can not be obtained analytically. They should instead be estimated using numerical procedures such as Monte Carlo (MC). Secondly, selecting the next point in EGO is a global maximization of the EI problem which has the same dimension as the design variables. In contrast, maximizing qEI and EEI has for dimension the number of the design variables multiplied by the number of new points (the number of free computing nodes). Therefore brute optimization of qEI and EEI using simple MC may not be cost effective. These problems have been partially addressed in [START_REF] Ginsbourger | Kriging is well-suited to parallelize optimization[END_REF], [START_REF] Ginsbourger | Dealing with asynchronicity in parallel Gaussian Process based global optimization[END_REF], [START_REF] Dubourg | Reliability-based design optimization using kriging surrogates and subset simulation[END_REF] by introducing heuristics. This report further studies expected improvement criteria for parallel and asynchronous computational environments. Firstly, the EI (µ,λ) criterion is formulated in a more general form. Secondly, in section 4, bounds on EI (µ,λ) are provided in section 3 and strategies for comparing EI (µ,λ) at various design points are discussed. Finally, the impact of the number of available computing units, λ, on the performance, i.e., the method speed-up, is investigated in section 5.

A unified presentation of parallel expected improvements

The traditional improvement in the objective value [START_REF] Donald | A taxonomy of global optimization methods based on response surfaces[END_REF], [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] at point x is defined as

I (ω) (x) = max(0, f min -min(Y (ω) )) = f min -min(Y (ω) ) + (1)
where f min = min(Y) is the minimum of the observations, x is the new point for the next simulation and Y (ω) = Y (ω) (x) |Y is the joint Gaussian vector (kriging metamodel) at the point of interest x conditioned on the past observations X, Y.

The multi-points improvement [START_REF] Ginsbourger | Kriging is well-suited to parallelize optimization[END_REF] can be rewritten as

I (λ) (ω) (x) = max(0, f min -min(Y λ (ω) )) = f min -min(Y λ (ω) ) + (2)
where f min = min(Y) is the minimum of the previous observations, x = (x 1 , ..., x λ ) are the λ new points for the next simulation and

Y λ (ω) = (Y 1 , ..., Y λ ) = Y (ω) (x 1 ), ..., Y (ω) (x λ
) |Y is joint Gaussian vector (kriging metamodel) at the new points x = (x 1 , ..., x λ ) conditioned on the past observations X, Y.

The multi-points asynchronous improvement is defined as

I (µ,λ) (ω) (x) = max(0, min(f min , Y µ (ω) ) -min(Y λ (ω) )) = min(f min , Y µ (ω) ) -min(Y λ (ω) ) + (3)
where f min = min(Y) is the minimum of the previous observations, x = (x µ+1 , ..., x µ+λ ) are the λ points for the next simulation and

Y λ (ω) = (Y µ+1 , ..., Y µ+λ ) = Y (ω) (x µ+1 ), ..., Y (ω) (x µ+λ ) |Y is again the Gaussian vec-
tor at the points of interest conditioned on the past observations X, Y, and

Y µ (ω) = (Y 1 , ..., Yµ) = Y (ω) (x 1 ), ..., Y (ω) (xµ)
|Y is the joint Gaussian vector at the points where the simulator has not yet provided responses (i.e., the busy points). The multi-points asynchronous improvement has first been introduced in [START_REF] Ginsbourger | Dealing with asynchronicity in parallel Gaussian Process based global optimization[END_REF], [START_REF] Ginsbourger | Dealing with asynchronicity in kriging-based parallel global optimization[END_REF]. As it can be readily seen from eq. (3), it is a natural measure of the progress made between points being or already calculated (the min(f min , Y µ (ω) ) term) and future points (the min(Y λ (ω) ) term). It has an important feature for asynchronous parallel calculation: it is null at the busy points. Indeed,

If min(Y µ (ω) ) ≤ f min , I (µ,λ) (ω) (x) = min(Y µ (ω) ) -min(Y µ (ω) ) + = 0 else min(Y µ (ω) ) > f min , I (µ,λ) (ω) (x) = f min -min(Y µ (ω) ) + = 0
When used in the context of optimization, a natural choice for the sampling criterion is to use the expectation of the improvement, thus for multi-points asynchronous improvement in eq. ( 3) we have

EI (µ,λ) (x) = E Ω (I (µ,λ) (ω) (x)). (4)
In the special case when µ = 0 and λ = 1, it reduces to classical EI, with a known analytical expression. When µ = 0 and λ ≥ 2, it is equivalent to the qEI and we will further write it as EI (λ) (note that for µ = 0 and λ = 2, an analytical expression is given in [START_REF] Ginsbourger | Kriging is well-suited to parallelize optimization[END_REF]). Note also that by applying the law of total expectations eq. ( 4) is equivalent to EEI [START_REF] Ginsbourger | Dealing with asynchronicity in parallel Gaussian Process based global optimization[END_REF]. In general calculating eq. ( 4) amounts to estimating a multidimensional integral (where the number of dimensions depends on the number of available and busy nodes) with respect to Gaussian density and must be based on numerical procedures, in particular MC methods.

MC based estimation procedure samples the Gaussian vector (Y λ , Y µ ) and estimates the criteria values. Even though the calculation of the improvement (eq. ( 2) and eq. ( 3)) is not complex, the estimation of the criteria may become time consuming. The error of MC estimate depends on the number of samples, which may be large if a high precision of the estimate is needed. Such cases are very likely to occur at the end of the optimization when the points x proposed by the mainly converged optimizer are close to each other. Additionally, the crude MC sampling becomes inefficient with the number of EGO iterations because the probability of improvement becomes small.

Illustration of criteria

For illustrating the criteria we use the one dimensional test function

y = sin(3x) -exp - (x + 0.1) 2 0.01 . ( 5
)
The test function together with predicted kriging mean and 2σ M C confidence intervals, one point EI, 2 points EI or EI (0,2) and 2 points asynchronous EI or EI (1,2) are shown in fig. 1. It can be well seen that the criteria are symmetrical with respect to the line

x 1 = x 2 . Also notice that the maximum of EI (0,2) (x) is close to x = (x 1 , x 2 ) where (x 1 , x 2 ) corresponds to the modes of EI(x).
The points obtained from max EI (λ) for several iterations are shown in fig. 2. In this case of 'deceptive' function, the region of the true optimum is located faster with increasing λ. Larger λ values induce a better exploration of the design space allowing construction of a globally more accurate regression model and thus preventing EI from stagnation. Figure 3 provides an example of points created through max EI (1,1) for 3 iterations (6 time steps), assuming that calculation of the objective takes two time steps, and that the second computing node becomes available after one time step.

EI (µ,λ) criteria bounds

In this section upper and lower bounds for the criteria of eq. ( 4) are derived. Firstly, in section 3.1, we address the parallel case EI (λ) and after, in section 3.2, we look at the asynchronous case EI (µ,λ) .

Bounds on the multipoints expected improvement

From the definition of multi-points improvement (c) Contour lines of 2 points EI,

I (λ) (ω) (x) = max(0, f min -min(Y 1 , ..., Y λ )) ≥ max(0, f min -Y i ) = I (λ=1) (ω) (x i ) (a)
EI (0,2) (x) (d) 3D plot of 2 points EI, EI (0,2) (x)
(e) Contour lines of 2 points async. EI (busy point x busy = -0.34), EI (1,2) (x) (f) 3D plot of 2 points async. EI (busy point

x busy = -0.34), EI (1,2) (x)
Fig. 1: Illustrations of the simple, the two-points and the two-points asynchronous Expected Improvements using the analytical test function eq. ( 5). EI (0,2) and EI (1,2) are calculated with 10000 MC simulations. Notice the maximum of EI (0,2) (x) is close to x = (x 1 , x 2 ) where (x 1 , x 2 ) corresponds to the modes of EI(x). Fig. 3: Example of points generated by max EI (1,1) during 6 time steps. It is assumed that calculation of the objective takes two time steps. Objective function (solid line), points used for kriging model (blue), points where the response value has been calculated (red), busy points (sent for calculation but not known yet -green). Initially and after the first time step θ = 0.3, for subsequent steps (after data arrives) θ is estimated by maximizing the kriging model likelihood.

and therefore the lower bound on multi-points EI is

EI (λ) (x) ≥ max i=1,λ EI (λ=1) (x i ). ( 6
)
The upper bound on the expectation of multi-points improvement can also be derived from its definition

EI (λ) (ω) (x) = E[max(0, f min -min(Y 1 , ..., Y λ ))] (7) = λ i=1 fmin -∞ +∞ yi . . . +∞ yi (f min -y i ) f N (µ,Σ) (y 1 , ..., y λ )(dy j ) j =i dy i ≤ λ i=1 fmin -∞ +∞ -∞ . . . +∞ -∞ (f min -y i ) f N (µ,Σ) (y 1 , ..., y λ )(dy j ) j =i dy i = λ i=1 fmin -∞ (f min -y i ) +∞ -∞ . . . +∞ -∞ f N (µ,Σ) (y 1 , ..., y λ )(dy j ) j =i dy i = λ i=1 fmin -∞ (f min -y i )f N (µi,σ 2 i ) (y i )dy i = λ i=1 EI (λ=1) (x i ) (8)
where f N (µ,Σ) is the density of a multivariate normal distribution with mean µ and covariance matrix Σ, (dy j ) j =i is the sequence from dy 1 to dy λ except dy i . The above relations provide the upper bound

EI (λ) (x) ≤ λ i EI (λ=1) (x i ). (9)
The plots in fig. 4 show upper and lower bounds versus EI (λ) (λ = 2) for 100 random points. The test case is the same as in the previous example of fig. 1 with the function eq. ( 5). EI (λ) is calculated using 100 and 10000 MC simulations. It can be seen that for these points bounds are seemingly tight.

Bounds on the asynchronous multi-points expected improvement

From the definition of asynchronous multi-points improvement, one has

I (µ,λ) (ω) (x) = max(0, min(f min , Y µ 1 , ..., Y µ µ ) -min(Y λ 1 , ..., Y λ λ )) (10)
≤ max(0, f min -min(Y λ 1 , ..., Y λ λ )) max(0, Y µ i -min(Y λ 1 , ..., Y λ λ )) (11)
and the expectation is

EI (µ,λ) (ω) (x) = ≤ E[max(0, f min -min(Y λ 1 , ..., Y λ λ ))] = EI (λ) (ω) (x) E[max(0, Y µ i -min(Y λ 1 , ..., Y λ λ ))] (12)
If we define and

I * (i,j) (ω) (x) = max(0, Y µ i -Y λ j ) (13)
EI * (i,j) (x) = E Ω (I * (i,j) (ω) (x)) (14)
then using a similar reasoning as in eq. ( 7) one can say that

E[max(0, Y µ i -min(Y λ 1 , ..., Y λ λ ))] ≤ λ j=1 EI * (i=1,j) (x j ). ( 15
)
The upper bound is

EI (µ,λ) (x) ≤ min   λ j EI (λ=1) (x j ), λ j EI * (i=1,j) (x j ), ..., λ j EI * (i=µ,j) (x j )   . (16)
Notice that all components on the right hand side of eq. ( 16) have analytical expressions since Y µ i Y λ j | Y is Gaussian with mean and variance known from kriging.

To estimate the lower bound on the asynchronous expected improvement one can write

I (µ,λ) (ω) (x) = max(0, min(f min , Y µ 1 , ..., Y µ µ ) -min(Y λ 1 , ..., Y λ λ )) = max(0, min(f min , Y µ 1 , ..., Y µ µ ) -Y λ 1 , ..., min(f min , Y µ 1 , ..., Y µ µ ) -Y λ λ ) ≥ (min(f min , Y µ 1 , ..., Y µ µ ) -Y λ j ) + .
The expectations of type E[(min(f min , Y µ 1 , ..., Y µ µ ) -Y λ j ) + ] are again integrals on the R µ+1 hyperspace truncated by several hyperplanes. As noted previously, the estimation of such integrals is based on numerical procedures and are not trivial. Instead a trivial lower bound of the asynchronous expected improvement from the definition of the improvement can be used

EI (µ,λ) (x) ≥ 0. ( 17
)
The plots in fig. 5 show upper bounds versus EI (µ,λ) (µ = 1, λ = 2) for a set of random points using the same test function eq. ( 5) as in the previous examples fig. [START_REF] Ginsbourger | Dealing with asynchronicity in parallel Gaussian Process based global optimization[END_REF] and fig. 1. EI (µ,λ) is calculated using 100 and 10000 MC simulations. It can be seen that there are clusters of points where the upper bound is considerably overestimated.

EI (µ,λ) estimation and selection of the candidate point

The optimization criteria of eq. ( 4) are intended to choose the most promising set of points for the next simulations. As already noted above, the computation of these criteria is not trivial, so that a simple heuristic strategy would be to use the bounds derived in the previous section. Because the bounds do not account for couplings between points (to the exception of the EI * values which account for pairwise couplings through Y µ i Y λ j in eq. ( 13)), the criteria values may be considerably under or over estimated (for multi-points EI (λ) fig. 4 and especially asynchronous multi-points EI (µ,λ) fig. 5). Therefore, we propose a mixed strategy which utilizes bounds together with MC estimations.

(a) Upper bound versus EI (µ,λ) for 100 random points with estimated 2σ confidence intervals. EI (µ,λ) estimates using 100 MC samples.

(b) Upper bound versus EI (µ,λ) for 1000 random points. EI (λ) estimates using 10000 MC samples.

(c) Lower bound (diamond), upper bound (square), estimated EI (µ,λ) (red dot) with estimated 2σ confidence intervals (red lines) on 100 random points. EI (µ,λ) estimates using 10 MC samples.

(d) Lower bound (diamond), upper bound (square), estimated EI (µ,λ) (red dot) with estimated 2σ confidence intervals (red lines) on 100 random points. EI (µ,λ) estimates using 100 MC samples. Fig. 5: Upper and lower bounds, EI (µ,λ) (µ = 1, λ = 2) for 100 and 1000 random points with 2σ confidence intervals. The test case is based on eq. ( 5) like in fig. 1. EI (µ,λ) is calculated using 100 and 10000 MC simulations. Notice in this example that for many points the upper bounds are overestimated.

Estimation and confidence intervals

The MC estimation allows us to estimate the expectations of improvements in eq. (4). By using crude MC the estimated EI (µ,λ) (denoted in this section as EI) is

EI M C (x) = 1 n N i=1 I (ω) (x)
where N is the number of samples and I (ω) (x) is the improvement calculated from the sampled trajectory of the conditional Gaussian process at the points of interest.

"Points" has a plural here because, as the reader may remember, in the case of the parallel asynchronous expected improvement, expectation is calculated over joint random variables located at µ + λ points in the optimization variable space. The error variance of this estimate also can be estimated from data as

σ 2 M C = 1 n(n -1) N i=1 (I (ω) (x) -EI M C ) 2 .
Therefore (EI M C |EI) follows Student's t-distribution with mean EI and variance σ 2

M C

with n -1 degrees of freedom. Bounds discussed in previous chapter provide us with prior information on EI. Lets assume that α ≤ EI ≤ β and EI uniformly distributed between its bounds, EI ∼ U (α, β). Using Bayes theorem it is possible to show that (EI|EI M C ) follows a truncated Student's t-distribution

p(EI|EI M C ) = p(EI M C |EI)p(EI) + -p(EI M C |EI)p(EI)dEI (18) = 1 σ M C fν EI -EI M C σ M C Fν α -EI M C σ M C -Fν β -EI M C σ M C -1 (19)
where Fν (.) and fν (.) are the c.d.f. and p.d.f. of the Student's t-variable with ν = n-1 degrees of freedom. The moments of truncated t-distribution are given for example in [START_REF] Kim | Moments of truncated student-t distribution[END_REF] and the formulas are copied in appendix A. Numerical tests have shown that these formulas become numerically unstable as the bounds are on the same side and far from the mean. However it is also know that as ν increases the t-distribution is well approximated by the Gaussian distribution. The density function of truncated Gaussian distribution is equivalent to eq. ( 18) except that instead of Fν (.) and fν (.) we have Φ(.) and φ(.) which are the c.d.f. and p.d.f. of the standard normal distribution. The mean and variance of truncated Gaussian random are known and if we denote by

u 1 = (α -EI M C )/σ M C and u 2 = (β -EI M C )/σ M C , then E(EI|EI M C ) = EI M C + φ(u 1 ) -φ(u 2 ) Φ(u 2 ) -Φ(u 1 ) σ M C (20)
and

VAR(EI|EI M C ) = σ 2 M C 1 + u 1 φ(u 1 ) -u 2 φ(u 2 ) Φ(u 2 ) -Φ(u 1 ) - φ(u 1 ) -φ(u 2 ) Φ(u 2 ) -Φ(u 1 ) 2 . ( 21
)
The

[E(EI|EI M C ) -k VAR(EI|EI M C ), E(EI|EI M C ) + k VAR(EI|EI M C )]
when k = 1 gives a 68% confidence that EI is inside this interval.

Ranking of candidate points

The confidence intervals from the previous section allow to compare two points. Let say that we have a set of points x i , i = 1..n , MC estimates EI M C (x i ) and σ 2 M C (x i ), let's denote by µ i = E(EI|EI M C (x i )) and σ i = VAR(EI|EI M C (x i )) . Now, with over 60% confidence (k = 1),

if µ i -kσ i ≥ µ j + kσ j then EI(x i ) ≥ EI(x j ) (22) if µ i + kσ i < µ j -kσ j then EI(x i ) < EI(x j ). ( 23
)
If neither eq. ( 22) nor eq. ( 23) hold then we do not have enough information to safely discriminate between the two points and additional data is necessary. In order to reduce the confidence intervals N = N i + N j additional MC samples are computed, where N i and N j are the number of samples at points x i and x j , respectively. A simple strategy is to add a number of samples proportional to the variance of the estimate,

N i = σ * i N σ * i + σ * j , N j = N -N i where σ * i = σ i , σ * j = σ j , if σ M C (x i ) > 0, σ M C (x j ) > 0.
Special care must be taken when none of the MC samples actually obtains trajectory where improvement is positive. In such cases the σ M C = 0 and EI M C = 0 are underestimated. On the other hand increasing number of total MC samples with no hits (trajectory providing positive improvement) indicates that true value of EI is small. To safely discriminate between two candidate points one needs to account for such situations.

If point x j has no MC hit but x i has at least one sample where improvement is positive it is possible to estimate the improvement variance σ 2 M C (x i ) and use this variance to overestimate the improvement variance at the other point x j as N i /N j σ 2 M C (x i ), where N i and N j are the number of already computed MC samples.Therefore if σ i > 0 and σ j = 0 then σ * i = σ i and σ * j is calculated using eq. ( 21) where σ M C (x j ) = σ M C (x i ) ni nj . However if no MC samples have positive improvement at both points .i.e. σ M C (x i ) = 0, σ M C (x j ) = 0 then the number of samples should be increased equally N i = N j = N/2, until the probability of improvement is small and we can assume that the true value of EI at both points is very small and further sampling is unnecessary, or until the MC budget is exceeded.

Illustration of the impact of bounds

In fig. 6 we show the 80th percentile of the number of MC simulations versus λ for comparing random pairs of points. For each λ setting, 100 random trajectories of Gaussian process are generated (and act as sample functions), the kriging models are built from 4 evenly spaced points using the same covariance form as that used in the process generation.

For each trajectory a pair of points is randomly selected (i.e., a set of µ points and two sets of λ points) and their EI (µ,λ) values are compared using MC based procedures with and without bounds. The initial number of MC samples is 10 and at each step N = 20. The maximum number of MC samples for each point is limited to 10 5 . If more MC samples are necessary one concludes that for the two points values of EI (µ,λ) are very similar and the points cannot be compared. The percentiles in fig. 6 are computed only on points that can be compared. For test cases where µ = 0, on the average 99 pairs of points where comparable within this number of maximum MC evaluations. For cases where µ = 1 and µ = 3, on the average only 75 and 69 pairs where comparable.

The results indicate that the bounds reduce the necessary number of MC evaluations, however (especially in asynchronous case µ > 0), the discrimination of points often requires a very large number of MC evaluations. Furthermore, in practice one is interested in selecting points with maximum EI (µ,λ) : the optimization of EI (µ,λ) eventually leads to regions where the values of EI (µ,λ) become similar and a much larger number of MC simulations are needed for the ranking of points i.e., the non comparable cases occur more often. Therefore the strategy based on mixing crude MC with bounds is effective only at the initial optimization steps when the global exploration of EI (µ,λ) resembles random sampling. Fig. 6: 80th percentile of the number of MC simulations necessary to discriminate between EI (µ,λ) at two random points versus λ. The number of MC simulations without bounds (dashed line), number of MC simulations with bounds (solid line).

5 An empirical study of EI (0,λ) scale-up properties To empirically compare the efficiency of the parallel improvement criteria (eq. ( 4)) to classical EI, we use random test functions in one dimension, i.e. we sample Gaussian process trajectories on a 1000 point grid using Matern 5/2 kernel functions [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF] with fixed scaling parameters θ = 0.3 (fig. 8, fig. 10) and θ = 0.15 (fig. 9, fig. 11 ).

Optimization strategy

For the maximization of EI (λ) (x), we use the Covariance Matrix Adaptation Evolution Strategy algorithm (CMA-ES, [START_REF] Hansen | The CMA evolution strategy: a comparing review[END_REF]) with some changes. Firstly, as we have data only at discrete locations, every point sampled by CMA-ES is mapped to the nearest grid point. Secondly, the ranking of points in a given population is based on the comparison procedure described in section 4. This MC based ranking makes the objective function noisy, which is compatible with the stochastic CMA-ES algorithm, but adds randomness in the tests.

When comparing two points x i and x j (two sets of points in actual space) which have close EI values (for example the points are close to each other) the necessary number of MC simulations for safe comparison may be very large. In order to keep the optimization procedure computationally feasible, the maximum number of MC simulation for estimating EI at one point is 10 5 . If the numbers of MC simulations for both points exceed this budget, we assume that it is not possible to discriminate between the two points (as explained in section 4) because their criteria values are very close. This situation occurs during the optimization, especially towards the end of the CMA-ES iterations, when the variance of the population is small. Furthermore, the values of EI (and probability of improvement) decreases with the number of iterations. This further increases the number of MC samples needed to discriminate between two points, as most of the sampled trajectories will be above the best observed point. These problems illustrate the necessity for more efficient criteria estimation procedures (such as MC variance reduction techniques), where more trajectories that are interesting for calculating the criteria ( below the best point) are used.

Results

The impact of λ on the speed-up of the parallel expected improvements is examined by comparing the actual improvements (difference between best objective values at the points of the initial DOE and at the λ points provided by criteria) after one time step.

The assumption is made that the simultaneous calculation of the objective is possible on multiple computing nodes, however the actual number of utilized processors depends on the criterion, i.e. EI (1) provides one new point, EI (2) two new points, etc.

The variation in actual improvement due to the EI optimization strategy is illustrated in fig. 7, where 10 runs of EI maximization are performed on a single test trajectory fig. 7a and the mean actual improvement with 1σ confidence bounds fig. 7c and mean improvement rank together with 1σ confidence fig. 7d. The improvement rank is obtained for each trajectory by ordering the best actual improvement values with respect to λ. It is used as a normalization method for better visualization of the speed-up related to parallel EIs. The mean improvement, mean rank with confidence intervals for actual improvement on 10 random test functions (fig. 8a and fig. 9a) at points of maxx EI (λ) (x) are illustrated in fig. 8c and fig. 9c. The initial kriging model is built using 4 evenly spaced points and the covariance hyper-parameter θ is fixed equal to the value used in trajectory generation. The results of this test indicate that the EI (λ) has a sub-linear improvement with respect to λ.

The λ impact on improvement is also studied in a more realistic scenario, where points of maximum EI are added sequentially to the DOE as they would during an optimization. The following typical three steps optimization strategy is investigated :

1. the kriging model is built from 4 initial points; 2. kriging covariance parameters are not known a priori and are estimated by maximizing likelihood; 3. EI (λ) (x) is maximized by CMA-ES; 4. the λ points together with f (x) are added to the DOE.

Steps 2, 3, 4 are repeated for 2 iterations.

The statistics of the overall best improvement (the best improvement after three iterations with respect to initial 4 point DOE) versus λ over 10 random test functions using the three step optimization strategy are shown in fig. 10 and fig. 11. As in the previous cases, the results suggest that the EI (λ) provides a sub-linear improvement with respect to λ.

Conclusion

In this study, the EI (µ,λ) criterion has been formulated by combining parallel and asynchronous versions of the expected improvement criterion. EI (µ,λ) measures the expected improvement brought by λ new points when the value of the objective function at µ already chosen points is not yet known but will be.

Bounds on EI (µ,λ) have been provided for use in MC based estimation and comparison strategy. The easily calculable bounds may be useful at the initial steps of EI (µ,λ) maximization when comparing points at distant locations and the coupling between design points is weak.

Finally, the effect of the number of available nodes, λ, on the optimization performance has been investigated in 1D test cases. The results indicate that parallel EIs provide sub-linear speed-up. The observed speed-up is in agreement with the intuition that the ability to obtain points in parallel should provide better and perhaps "safer" results. Note that in the test cases studied here the actual type of kriging covariance kernel is known. It seems, however, that in general the parallel calculation of λ points provides better exploration of the design space. Therefore it may partially solve the problem of deceptive functions [START_REF] Donald | A taxonomy of global optimization methods based on response surfaces[END_REF], [START_REF] Forrester | Global optimization of deceptive functions with sparse sampling[END_REF]. Better exploration of the design space allows globally more accurate construction of kriging models and therefore reduces the risk of stagnating searches in local regions due to deceptive initial states.

The implementation of the EI (µ,λ) (x) for optimization indicates several difficulties. Firstly, the necessary number of MC simulations increases a lot in order to be able to discriminate close points. Secondly, crude MC sampling is often inefficient because trajectories that are better than the best observation become seldom as the optimization proceeds. Thirdly, the number of dimensions of x increases proportionally to the number of nodes. These challenges call for further studies on specialized optimization heuristics and MC procedures.

  Test function eq. (5) (dashed), kriging model based on 5 observations -predicted mean (blue), 2σ predicted confidence intervals (red). (b) One point EI, EI(x).

Fig. 2 :

 2 Fig. 2: Example of points generated by max EI (λ) (red) at different iterations, true function (solid line), points used for kriging model (blue). First line λ = 1, second line λ = 2, third line λ = 3, forth line λ = 4. First column -first iteration, second column -second iteration, third column -third iteration. For the first iteration θ = 0.3, for iterations 2 and 3, θ is estimated by maximizing the kriging model likelihood.

  (a) Upper bound versus EI (λ) with estimated 2σ confidence intervals. EI (λ) estimated using 100 MC samples. (b) Upper bound versus EI (λ) with estimated 2σ confidence intervals. EI (λ) estimated using 10000 MC samples. (c) Lower bound versus EI (λ) with estimated 2σ confidence intervals. EI (λ) estimated using 100 MC samples. (d) Lower bound versus EI (λ) with estimated 2σ confidence intervals. EI (λ) estimated using 10000 MC samples.

Fig. 4 :

 4 Fig.4: Upper and lower bounds versus EI (λ) (λ = 2) for random chosen points with 2σ confidence intervals. The test case is based on eq. (5) as in fig.1. EI (λ) is calculated using 100 and 10000 MC simulations. Notice that for the points in this example bounds are seemingly tight.

  (a) µ = 0. Sample trajectories generated on 1000 point grid with Matern 5/2 kernel and θ = 0.3. (b) µ = 0. Sample trajectories generated on 1000 point grid with Matern 5/2 kernel and θ = 0.15. (c) µ = 1. Sample trajectories generated on 1000 point grid with Matern 5/2 kernel and θ = 0.3. (d) µ = 3. Sample trajectories generated on 1000 point grid with Matern 5/2 kernel and θ = 0.3.

  (a) Sample trajectory on 1000 point grid with Matern 5/2 kernel and θ = 0.3. (b) Actual improvement of EI (λ) for 10 independent maximizations, λ = 1...4. (c) Actual mean improvement, 1σ confidence intervals of EI (λ) , λ = 1...4. (d) Mean improvement rank, 1σ confidence intervals of EI (λ) , λ = 1...4.

Fig. 7 :

 7 Fig.7: Variation in actual improvement due to the stochastic maximization of EI(λ) . Statistics calculated for 10 runs on a fixed test function.

  (a) Sample trajectories on 1000 point grid with Matern 5/2 kernel and θ = 0.3. (b) Actual mean improvement over 10 trajectories for λ = 1...4. (c) Mean improvement rank, 1σ confidence intervals over 10 trajectories for λ = 1...4.

Fig. 8 :

 8 Fig. 8: λ impact on performance observed on 10 relatively smooth (θ = 0.3) random functions

Fig. 9 :

 9 Fig. 9: λ impact on performance observed on 10 shaky (θ = 0.15) random functions

Fig. 10 :

 10 Fig. 10: λ impact on performance observed on 10 relatively smooth (θ = 0.3) random functions over 3 iterations.

  (a) Sample trajectories on 1000 point grid with Matern 5/2 kernel and θ = 0.15. (b) Actual mean overall improvement after 3 iterations, over 10 random trajectories, λ = 1...4. (c) Mean overall improvement rank, 1σ confidence intervals over 10 random trajectories, λ = 1...4.

Fig. 11 :

 11 Fig. 11: λ impact on performance observed on 10 shaky (θ = 0.15) random functions over 3 iterations.

A Moments of Students t-truncated distribution

If we denote by

] are given for example in [START_REF] Kim | Moments of truncated student-t distribution[END_REF] and

-bB

)

where


Yann Gavet

Mathieu Fernandes

Jean-Charles Pinoli

Automatic quantitative evaluation of image registration techniques with the dissimilarity criterion in the case of retinal images

Keywords: Image registration, Human retina, Dissimilarity criterion

In human retina observation (with non mydriatic optical microscopes), a registration process is often employed to enlarge the field of view. For the ophthalmologist, this is a way to spare time browsing all the images. A lot of techniques have been proposed to perform this registration process, and indeed, its good evaluation is a question that can be raised. This article presents the use of the dissimilarity criterion to evaluate and compare some classical featurebased image registration techniques. The problem of retina images registration is employed as an example, but it could also be used in other applications. The images are first segmented and these segmentations are registered. The good quality of this registration is evaluated with the dissimilarity criterion for 25 pairs of images with a manual selection of control points. This study can be useful in order to choose the type of registration method and to evaluate the results of a new one.

INTRODUCTION

An image registration technique consists in representing two or more images into the same coordinate system, such that there is a correspondance between the structures of the images. For more informations, a survey is presented in [START_REF] Zitova | Image registration methods: a survey[END_REF]. When performing an image registration, the question of accuracy is raised. A visual assessment can be performed, but there is not quantification. Some similarity measures can be used [START_REF] Ritter | Registration of stereo and temporal images of the retina[END_REF][START_REF] Skerl | A protocol for evaluation of similarity measures for rigid registration[END_REF][START_REF] Lalibert | Registration and fusion of retinal images-an evaluation study[END_REF]. This article proposes the use of a novel dissimilarity criterion [START_REF] Gavet | The criterion: a contour-based supervised tolerant image segmentation evaluation[END_REF] to evaluate quantitatively the performance of different registration methods. This article is centered on this dissimilarity criterion and is not intended to be a survey on registration techniques. The problem of the quality control on retinal images registration is addressed as an illustration, but the proposed dissimilarity criterion can be employed in a lot of different application areas such as the the control of image segmentation processes [START_REF] Gavet | The criterion: a contour-based supervised tolerant image segmentation evaluation[END_REF].

Registration techniques can be divided into two categories: iconic registration and feature-based registration techniques. In the case of retinas, iconic registration is not a good choice because images present different illumination conditions. Moreover, one can notice that the ophthalmologists manually perform this registration process by using a landmark matching based on the vasculature and its bifurcations 6.

Dissimilarity criterion

The dissimilarity criterion is defined as a comparison either between continuous [START_REF] Gavet | The criterion: a contour-based supervised tolerant image segmentation evaluation[END_REF] sets or between discrete [START_REF] Gavet | The criterion: a contour-based supervised tolerant image segmentation evaluation[END_REF] sets (in this article, the considered sets are discrete, e.g. finite sets of pixels). Practically, these sets can be the results of some image segmentation process. This criterion introduces a tolerance in the measure (thus, it is not mathematically a metric) by the mean of the Minkowski addition.

Further author information: (Send correspondence to Y. Gavet) Y. Gavet: E-mail: gavet@emse.fr, Telephone: [33] 4 77 42 01 70

Minkowski addition

The Minkowski addition [START_REF] Minkowski | Volumen und Oberfläche[END_REF] defines an algebraic operation between sets in R n .

If M and B are two sets of R n , the Minkowski sum of M and B, denoted M ⊕ B, is then defined by:

M ⊕ B = {m + b|m ∈ M, b ∈ B} M ⊕ B = b∈B {m + b|m ∈ M } ( 1
)
where ⊕ is the Minkowski addition symbol. In the field of Mathematical Morphology [START_REF] Soille | Morphological Image Analysis: Principles and Applications[END_REF], it corresponds to the dilation, and B is called a structuring element (for example a ball).

Definition of the dissimilarity criterion

The dissimilarity criterion with the tolerance ρ applied to segmented (binary) images is defined in the case of discrete (numerical) images by the following equation [START_REF] Gavet | The criterion: a contour-based supervised tolerant image segmentation evaluation[END_REF]:

ρ M (X) = #{(X \ M ⊕ ρN ) ∪ (M \ X ⊕ ρN )} #{M ⊕ ρN } (2)
with N being the structuring element of radius 1 (the unit ball) and # denoting the number of pixels in the set (# is the cardinal operator, counting the number of non null pixels in the sets X or M . The reader can see this formula as a tolerant version of the symmetric difference). Practically, ρ is the radius of the ball used to dilate the binary images. The choice of this value is not the purpose of this article, though it remains an important question. The value ρ = 3 is chosen.

METHOD

Human retina image database

A database of 25 pairs of retina images has been constructed. For each pair of images, the points of correspondance (called control points) have been manually marked by an expert (a minimum of 15 points of correspondance for each pair). Generally, these points are automatically extracted [START_REF] Fernandes | A feature-based dense local registration of pairs of retinal images[END_REF], but this remains a difficult task. The manual selection of the pairs of points thus avoids the difficult problem of their automatic detection.

Additionnaly, the (almost circular) region of interest of the images are extracted (see Fig. 1): this will be used as a mask to constrain the overlapping areas.

Image segmentation

A segmentation process extracts the vessels present in the images. A lot of different methods exists for achieving this task [START_REF] Chaudhuri | Detection of blood vessels in retinal images using two-dimensional matched filters[END_REF][START_REF] Staal | Ridge-based vessel segmentation in color images of the retina[END_REF][START_REF] Chanwimaluang | Hybrid retinal image registration[END_REF]. In this article, only one method is employed to get the vessels: the one of Chanwimaluang [START_REF] Chanwimaluang | Hybrid retinal image registration[END_REF]. Two examples are presented in Fig. 1(c) and 1(d). These segmented vessels will be used later for the quantitative evaluation of some registration techniques.

As we choose to evaluate the registration process by the evaluation of the registration of the segmentation, it is important to notice that some other segmentation evaluation criterions could have also be used 5. The main drawbacks of these methods is that they are not tolerant to small spatial variations of the segmentated objects.

Image registration

The registration process consists in finding a spatial transformation to be applied to one image to fit the other. Thus, the computation of this transformation is performed on the sets of correspondance points. If u and v are the coordinates of a point in the image to be registered U (as Unregistered), and x and y its coordinates in the second image B (as Base image), the transformation used belongs to one of the following categories presented hereafter (see [START_REF] Zitova | Image registration methods: a survey[END_REF][START_REF] Imgfsr | Transformation functions for image registration[END_REF] for more details on these methods): • Similarity: the computed transformation is a similarity (a combination of translation (h and k), rotation θ and scaling s).

u = x • s cos θ -y • s sin θ + h v = x • s sin θ + y • s cos θ + k
• Affine: a first order polynomial relation between the two systems of coordinates.

u = ax + by + c (3) v = dx + ey + f (4)
• Projective: this transformation may be used in the case of nonlinearities in the lens or sensor.

u = ax+by+c dx+ey+1 (5) v = fx+gy+h dx+ey+1 ( 6
)
• Polynomial: the relation between the coordinates is polynomial. Let P and Q be two polynomials for which the maximum degree in x and y is n. Practically, n is 2 or 3. If n = 1, the transformation is affine.

u = P (x, y) ( 7 ) v = Q(x, y) ( 8 )
• Piecewise linear [START_REF] Goshtasby | Piecewise linear mapping functions for image registration[END_REF]: a linear transformation function is used to map a triangle (defined by a Delaunay triangulation of the control points) in the registered image to the corresponding triangle in the base image.

• Local linear: Fernandes et al. [START_REF] Fernandes | A feature-based dense local registration of pairs of retinal images[END_REF] proposed a method that computes a transformation for each point in the images from its two nearest control points.

The masks of the images are applied so that only the overlapping regions appear as a result.

QUANTITATIVE EVALUATION AND RESULTS

The goal of this article, presented in this section, is to propose an answer to the questions: can I control the good quality of a registration process, how good is it, and how can I choose among the different transformations?

Mean results over the database

The registration process is performed on the segmented images of human retinas. Only the overlapping regions are kept (with the use of the masks) and the dissimilarity criterion is applied on them. To compare the different methods together, the results are normalized by the results of the affine method, which is considered as a reference. Notice that the lower the criterion value is, the better the registration method is.

Let B be the base (segmented) image, and U the unregistered (segmented) image. The transformation T maps U into the coordinate system of B. If A is the affine transformation, the normalized criterion, denoted ˜ , is computed as follows:

˜ B (M ) = 3 B (T (U )) 3 B (A(U )) (9)
• ˜ ≤ 1 means that the transformation gives a better result than the affine transformation.

• ˜ ≥ 1 means that the transformation gives a lower result than the affine transformation.

The results are proposed in Fig. 3 and Tabs. 2 and 1. They show that the second order polynomial method appears better than the others. The local linear method would be a good choice, but present high values in some cases.

These results must be tempered by the fact that they depend on the choice of the control points, this is why the next subsection will try to evaluate their influence.

Registration method

Influence of the number of control points

To evaluate the influence of the number of control points, one pair of corresponding human retina images is taken (Pair number 25). This pair of images contains 22 pairs of control points. The question we want to answer is: what happens if we take only a few control points among those?

The results are presented in Fig. 4. The number of control points is randomly choosen (among 22 points) and the average value is represented is abscissa 100 random choices. 3 B (T (U )) (measuring the dissimilarity of the transformation of the unregistered image T (U ) compared to the base image B) is represented in ordinates. Notice that each method requires a minimum number of pairs of control points.

It appears (for this pair of images) that the local linear method gives better results even for a low number of control points. This is only when using more than 17 pairs of points that the second order polynomial method gives lowest and thus better results. . This figure shows the evaluation of the registration in regard to the number of control points. This example is computed for only one pair of retina images. It shows that the so-called local linear method generally performs well even when a small number of points is used.

CONCLUSION AND PERSPECTIVES

The problem of retinal images registration is adressed as an illustration to the use of the dissimilarity criterion. In this article, some classical registration techniques are evaluated and compared. It is shown that the criterion is fairly simple to understand and manipulate. The reader can notice that could also be employed to evaluate results of segmentation, which was the reason why it was originally introduced [START_REF] Gavet | The criterion: a contour-based supervised tolerant image segmentation evaluation[END_REF].

Figure 1 .

 1 Figure 1. Example of a pair of retina images. The represented images (a) and (b) are the green component of the RGB images. The binary result of the segmentation (see subsec. 2.2) of these images is respectively represented in (c) and (d) for the images (a) and (b).

Figure 2 .Figure 3 .

 23 Figure 2. Result of the registration of segmented images of a retina. The masks of the regions of interest have been applied. The base image is represented by the red channel, the registered image is represented by the green channel. The superimposition appears in yellow.

3 B

 3 (T (U ))

Figure 4

 4 Figure 4. This figure shows the evaluation of the registration in regard to the number of control points. This example is computed for only one pair of retina images. It shows that the so-called local linear method generally performs well even when a small number of points is used.

Table 1 .

 1 Results over the 25 images of the retina database. The best average value is obtained by the second order polynomial method. The local linear method gives generally good results but sometimes the distortion of the registered images is too strong (see Fig.3).

		Mean ˜
	Affine	1
	Similarity	0.99
	Projective	0.93
	Polynomial of order 2	0.68
	Polynomial of order 3	0.84
	Piecewise linear	0.95
	Local linear	0.86

Table 2 .

 2 This table presents the normalized results (˜ ) for the 25 pairs of retina images and for the different registration methods. The lower values are the best results. The first column presents the number of pairs of control points for each pair of retina images in the database.

ACKNOWLEDGMENTS

The authors wish to thank Pr. Gilles Thuret and Pr. Philippe Gain from University Hospital Centre, Saint-Etienne, France for supporting this work and for providing pictures of human retinas.


Ola Suleiman Ahmad
email: ahmad@emse.fr

Johan Debayle

Nesrine Gherras

Benoît Presles

Jean-Charles Pinoli

Gilles Févotte

Ola Suleiman

Recognizing overlapped particles during a crystallization process from in situ video images for measuring their size distributions

Keywords: Crystallization, Salient Corner Detection, Grouping Method, Shape Recognition, Particle Size Distributions

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Crystallization is a chemical engineering process which is used in time-continuous mode in order to produce solid products such as sodium, ammonium sulfate and potassium chlorides, and in batch mode to fabricate fine chemicals and drugs. It consists of two major events. The first event is the nucleation, where the solute molecules distributed in the solvent aggregate and increase the concentration of the solute, then stabilize under certain conditions to form the nuclei. The second event is the crystal growth, where the nuclei growths and reaches a critical size to finally form one or different crystal structures called polymorphisms. Controlling the particle shape and size distributions (PSD) is very important to precisely evaluate the quality of the final drugs and chemical products, and to the downstream processes. [START_REF] Anda | Classifying organic crystals via in-process image analysis and the use of monitoring charts to follow polymorphic and morphological changes[END_REF] The main challenge to achieve a successful and robust characterization of the PSD by image analysis is the segmentation process, which aims to separate the particles from the background and ones from each others. It becomes a critical problem in crystallization because of several reasons, such as: the perspective projection of the 3D crystal shape onto the image plane, the blurred appearance of unfocused particles, the heterogeneity and the anisotropy of the solute, the degree of agglomeration or overlap, and the random variation in size of the observed particles. Many methods have been proposed to detect the crystals for monitoring the particle size distribution. An algorithm 2 based on line segment detection 3 was proposed to detect rectangular crystals by assuming that the rectangular object can be constructed from two or more line segments. In their preceding work, [START_REF] Larsen | Model-based object recognition to measure crystal size and shape distributions from in situ video images[END_REF] the algorithm is modified to be used to detect prismatic crystals basing on 3D-model fitting. In both applications, the results was satisfied in detecting accurately the overlapped particles from images but in low solids concentrations, where the particles are small enough with no high degree of overlap. However, the performance of the algorithm declines in high solids concentrations and with the existence of unfocused objects. A multi-scale segmentation method 5 was proposed and found to be effective for extracting crystals from unclear background of varied illumination, but the overlapped particles are detected as one object which is unwanted for the present paper setting. In the recent work 6 related to the last one, 5 the concave points are extracted from the detected particles 5 which are defined as markers for a marker-controlled watershed algorithm 7 in order to separate the overlapped crystals.

The method showed efficiency in separating the attached particles even at high solids concentrations, but it was not sufficient to detect the overlapped crystals of polygonal shapes. A restoration method 8 was proposed for detecting unfocused crystals from a segmented image [START_REF] Beucher | The watershed transformation applied to image segmentation[END_REF] for achieving an accurate measurement of the particle size distribution. However, the method detects the overlapped crystals as individual ones.

The goal of this paper is to outline an applied efficient method, based on image analysis, 9, 10 that will be able to separate the overlapped particles of rectangular and regular/irregular shapes during the batch crystallization process at different solids concentrations. This paper is organized as follows. The experimental setup of the crystallization process is presented in section 2. The image segmentation method of overlapped particles is introduced in section 3. The application of this method for analyzing the particle shape and size distributions is described in section 4. Finally, a conclusion and perspectives are drawn in section 5.

EXPERIMENTAL SETUP

The experimental setup of the batch crystallization process of Ammonium Oxalate crystals is conducted in pure water. [START_REF] Gherras | Observations sur les nucléations multiples en présence d'impuretés et conséquences pratiques[END_REF] Figure 1(a) illustrates the materials and the imaging system of the experiment. The reactor is a 2L, double-jacketed, glass vessel provided by a condenser. The mixer is a 3 blades, profiled, propeller, and baffles (Mixell TT) driven by a motor with a speed of 300 rpm which is sufficient to obtain a homogeneous suspension. The temperature of the crystallizer is controlled by a heating/cooling bath provided with an external probe of temperature P t100Ω in a feedback control way. The In situ imaging system consists of a probe (EZProbe sensor) [START_REF] Presles | Novel image analysis method for in situ monitoring the particle size distribution of batch crystallization processes[END_REF][START_REF] Presles | Monitoring the particle size distribution using image analysis during batch crystallization processes[END_REF] equipped with a CCD monochrome camera with high resolution up to 4μm 2 per square pixel, and 640 × 480 pixels per image. The video images are acquired by a frame grabber with a rate of 25 fps. Images of particles are formed using a transmitted light conducted within the probe by optical fibers from the light source. Figure 1(b) illustrates an in situ image example of Ammonium Oxalate crystals obtained during this experiment.

IMAGE SEGMENTATION METHOD

A method based on the geometric properties of the particle shapes has been developed to detect the overlapped rectangular and regular/irregular prismatic crystals by the assumption that each particle can be constructed from grouping three of their correspondent salient corners that belong to the same particle. A computer simulation process has been presented in order to validate the method, and then experiments have been implemented on the real images. The method (see Fig. 2) is divided into two main parts involving the detection of salient corners from the gray scale image, then recognizing the shape by clustering the correspondent salient corners that satisfy certain geometric conditions.

Salient Corner Detection

Many approaches have been developed for the corner detection because of its importance in computer vision, and pattern recognition applications. Some methods [START_REF] Harris | A combined corner and edge detection[END_REF][START_REF] Smith | A new approach to low level image processing[END_REF][START_REF] Rosten | Fusing points and lines for high performance tracking[END_REF] extract the corners directly from the gray level image. However, they are not sufficient for the present work since the orientation of salient corners is also required (see Eq. 1). Other methods [START_REF] Zhang | Robust image corner detection based on scale evolution difference of planar curves[END_REF][START_REF] Shen | Corner detection based on modified hough transform[END_REF][START_REF] Teh | On the detection of dominant points on digital curves[END_REF] extract the boundaries of the objects and then detect the corners from the planar curves or arc segments. Although, the performance of some of these methods depends on the robustness of extracting the arc segments from the gray level image. [START_REF] Shen | Corner detection based on modified hough transform[END_REF] A salient corner C i is defined, in this work, by the intersection of its two associated edges in the image spatial support D as follows:

C i = (x, y, u 1 , u 2 ) i : (x, y) ∈ D, u 1 2 = 1, u 2 2 = 1 (1)
The salient corner detection method is based on extracting the candidate corners from the gray level image, and then on analyzing these corners in the edge image (see Fig. 2) to finally identify the salient corners. Firstly, a preprocessing on the gray level image is needed to increase the contrast and to remove the illumination heterogeneity of the background (see Fig. 3(c)), then the candidate corners are detected from the restored image by Rosten method [START_REF] Rosten | Fusing points and lines for high performance tracking[END_REF] as illustrated in Figure 3(d). Secondly, the binary edge image is obtained using a multi-scale edge detection, [START_REF] Anda | Multi-scale segmentation image analysis for the in-process monitoring of particle shape with batch crystallisers[END_REF] because it shows efficiency in extracting the boundaries of the semi-transparent overlapped particles whereas the strength of the edge is affected by the overlap (see Fig. 3(e)). Canny detector [START_REF] Canny | A computational approach to edge detection[END_REF] is used to get the edges from different scales, then a set of morphological operations [START_REF] Soille | Morphological image analysis: principles and applications[END_REF] and edge cleaning are done to extract the final contour of the overlapped particles (see Fig. 3(f)). Finally, the direction of the salient corners is obtained from the directional correlation function, denoted E, which is defined in Equation 2:

∀θ ∈ [0, 5 • , 10 • , • • • , 360 • [ : E (θ) = M x=1 N y=1 W Ci (x, y) × T θ (x, y) ( 2 )
where W Ci is the binary edge image restricted to a window of size M × N pixels centered at corner C i , and T θ is a line template oriented at θ • . The correlation function of the salient corners has two major peaks that refer to its associated edges. Other corners such as junctions will have more than two dominant peaks, while the correlation values at the isolated points will be near zero with no significative peaks. Consequently, all corners which have less or more than two dominant peaks are removed. So an image of salient corners is obtained as illustrated in Fig. 3(g)

Shape Recognition from clustering method

Following the salient corner detection, a clustering technique is performed for extracting the final shape of the particles. The procedure depends on searching three correspondent corners that satisfy geometric conditions and thus belong to the same particle. The clustering technique is developed for recognizing rectangular and regular/irregular prismatic overlapped particles, whereas the general shape of the Ammonium Oxalate crystals is shown to be like rectangular or prismatic shape with faces of different growth rate under the experimental conditions. [START_REF] Jakubczyk | Morphology and kinetic of crystal growth[END_REF] The method is evaluated on 150 images acquired in low, intermediate and solids concentrations, successively (see Fig. 4) and a quantitative comparison with the visual recognition by a human expert is illustrated in Table 1.

The results show that the method is efficient in detecting the particles in low, intermediate, and even high solids concentrations. Indeed, more than 75% of the particles could be sightly recognized, which is sufficient for monitoring the size distributions.

PARTICLE SIZE DISTRIBUTIONS

Following the segmentation process, the particle size distributions are illustrated for 3 sets of video sequence images, each set contains 50 images, which are acquired in low, intermediate and high solids concentrations. The length and width of each detected crystal is calculated for all image sets and the size distributions are evaluated [START_REF] Scott | On optimal and data-based histograms[END_REF] for illustrating the frequency of each length and width parameters of the crystals during the experiment (see Fig. 5). The results show that the particles have a tendency to exhibit a prismatic shape in high agglomeration, whereas the small particles exhibit a rectangular shape. Small particles refer to the nucleation process and the beginning of the experiment. Thus, the Ammonium Oxalate crystals have a rectangular shape in low concentrations, then the faces grow to take a prismatic shape in high concentrations. This analysis is important for characterizing the flow properties, dissolution rate, agglomeration tendencies, mixing and demixing properties, filterability, etc. [START_REF] Allen | Powder sampling and particle size determination[END_REF] A quantitative comparison between automatic and manual sizing is performed (see Fig. 6) using the cumulative distribution function (CDF), which is independent of the bins width or number contrary to the particle size distribution function. The normalized error is calculated for both length and width cumulative distributions: where Cumh M , Cumh A are the cumulative distributions of manual and automatic sizing, respectively. The results are satisfying and show that the proposed automatic recognition method could be used for analyzing the particle size distributions of the overlapped particles.

Err = |Cumh M -Cumh A | Cumh M ( 3)

CONCLUSION AND PERSPECTIVES

A recognition method of overlapped particles based on salient corner detection and geometric grouping has been described and applied. The method is proposed to recognize rectangular, regular and irregular prismatic shapes. Then, the particle size distributions are measured on a set of video sequence images in different solid concentrations. A quantitative comparison between the present method and the visual recognition by a human expert is performed. The results show that the method can recognize generally more than 96%, 81% and 75% of the particles recognized by the human vision in low, intermediate and high solids concentrations, respectively. The cumulative distribution functions of particles width and length are calculated for both manual and automatic sizing. The maximum normalized difference between the cumulative distribution functions was 1.5%. The method can be generalized to detect parallelograms or other convex quadrilateral forms. It could be modified for recognizing other geometric polygons by clustering their correspondent salient corners according to the geometric properties of the desired shape. In the difficult cases of overlap, when the image contains a large number of semitransparent overlapped particles, it becomes difficult to recognize accurately the particles, which might decrease the recognition percentage. So, our future work is to ameliorate the method by estimating the shape using the likelihood function and the clustering propability of two or three salient corners that belong to the same polygonal shape.

Figure 1 .

 1 Figure 1. (a) Experimental setup of batch crystallization process. (b) A gray level image of Ammonium Oxalate crystals in pure water.

Figure 2 .

 2 Figure 2. The diagram of the proposed method

Figure 3 .Figure 4 .

 34 Figure 3. Example of the salient corner detection applied on a real image of Ammonium Oxalate crystals. (a) Original gray level image. (b) Background image. (c) Removing the background from image (a). (d) Candidate corners image. (e) Multi-scale edge detection. (f) Edge image after contour completion and edge cleaning steps. (g) Salient corners image.

Figure 5 .

 5 Figure 5. Particle length and width distributions. (a) Particle width distribution of prismatic and rectangular crystals. (b). Particle length distribution of prismatic and rectangular crystals.

Figure 6 .

 6 Figure 6. Comparison of cumulative distribution functions (CDF) of particles length and width for 150 in-situ images. (a) Comparison of particles width distributions. (b). Comparison of particles length distributions.

Table 1

 1 Quantitative comparison between automatic and visual recognition for 150 images of overlapped Ammonium Oxalate particles:

	Particles (numbers)	Low overlap Intermediate overlap High overlap
	Manually recognized (N V )	130	149	374
	Automatically recognized (N A )	125	137	327
	Right N R	125	122	281
	Missed N M = N V -N A	5	12	47
	False N F = N A -N R	0	15	46
	Right to automatic recognition ratio
	in percentage (N R /N A )%	100%	89.1%	85.9%
	Missed to visual recognition ratio
	in percentage (N M /N V )%	3.9%	8.1%	12.6%
	False to automatic recognition ratio
	in percentage (N F /N A )%	0%	10.9%	14.1%
	Right to visual recognition ratio
	in percentage (N R /N V )%	96.2%	81.9%	75.2%

ACKNOWLEDGMENTS

Financial support from French Research Agency ANR is partly provided for the project IPAPI (Improving the Properties of Active Pharmaceutical Ingredients), ref. 07-BLAN-0183.


Sylvie Chambon
email: chambon@ifsttar.fr

Alain Crouzil
email: crouzil@irit.fr

COMBINATION OF CORRELATION MEASURES FOR DENSE STEREO MATCHING

Keywords: Stereovision, matching, correlation, classic measures, robust statistics, fusion

In the context of dense stereo matching of pixels, we study the combination of different correlation measures. Considering the previous work about correlation measures, we use some measures that are the most significant in five kinds of measures based on: cross-correlation, classic statistics, image derivatives, nonparametric statistics and robust statistics. More precisely, this study validates the possible improvement of stereo-matching by combining complementary correlation measures and it also highlights the two measures that can be combined in order to take advantage of the different methods: Gradient Correlation measure (GC) and Smooth Median Absolute Deviation measure (SMAD). Finally, we introduce an algorithm of fusion that allows to combine automatically correlation measures.

INTRODUCTION

Finding homologous pixels in a stereo pair of images is one of the most important step in order to recover the 3D structure of a scene by stereovision. Many methods have been proposed in the literature where local methods are distinguished from global ones. More precisely, matching methods can be described with essential components, this term has been firstly introduced by [START_REF] Scharstein | A taxomomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF]. These components are: the matching cost, the optimization method, the introduction of multiple passes, i.e. to improve the matching performances, some approaches are based on several methods applied in sequence. This description leads to a four type categorization: local methods, global ones (without correlation measure), mixed method (global ones with a correlation measure) and the methods with multiple passes. Our purpose is to introduce a multipass algorithm based on combination of local methods.

Local methods are easy to implement, low time consuming, quite efficient, and consequently intensively used.

Unfortunately, characterising how the existing correlation measures are effective, i.e.

obtaining correct matching in different areas of the image, is still an open issue.

In our work [START_REF] Chambon | Occlusions handling in dense stereo matching[END_REF] on local costs, the influence of different measures on the quality of stereo matching results have been studied, in particular, near occluded regions and, in [START_REF] Chambon | Towards correlationbased matching algorithms that are robust near occlusions[END_REF], we demonstrated that a measure based on a robust statistics tool combined with a cross correlation measure allows to obtain better performances than using a correlation measure alone. These results raise up three new questions:

• Which are the correlation measures that are the most complementary to cover all the matching difficulties?

• Is it advantageous to combine numerous measures and how many?

• Following up the previous questions, can we propose an algorithm that combines more than one measure to obtain a whole dense and correct matching (or disparity) map1 and is it more efficient than the method based on a sole measure?

Existing methods are briefly presented before the description of the data set used for validating the proposed method. Then, combination study is described leading to the proposal for matching algorithm based on merging the results obtained from various correlation measures. Finally, results are presented.

CORRELATION MEASURES

The principle of a local cost, i.e. a correlation measure, is to consider that two homologous pixels and their respective neighborhoods, are similar, from a photometric point of view. The main difficulties of these methods are: illumination changes, untextured areas and occlusions. Many measures have been introduced to tackle out these difficulties. Based on the results of 41 measures on a benchmark of 42 images, presented in [START_REF] Chambon | Occlusions handling in dense stereo matching[END_REF], we propose to study the complementarity of these measures, and, in particular, the best measures of each families.

I w

The images with w ∈ {l, r} (left and right).

I i, j w , p i, j w
The grey level of the pixel p i, j w of coordinates

(i, j) in image I w is I i, j w . Moreover, p x,y r is the correspondant pixel of p i, j l . N *
The number of pixels in the neighborhood:

N f = (2N v + 1) ×(2N h + 1), N v , N h ∈ N * . f w
The vector of grey levels of pixels in the correlation windows (in I w ):

f w = (• • • I i+p, j+q w • • • ) T = (• • • f k w • • • ) T where T is the matrix transposition operator and p ∈ [-N v ; N v ], q ∈ [-N h ; N h ].
f w The mean of the grey levels in f w . f k w The element k of vector f w .

L P

The L P norms:

f w P = (∑ N f -1 k=0 | f k w | P ) 1/P with P ∈ N * and f w = f w 2 .
In the following description, when no explicit reference is given, the reader should consult [START_REF] Aschwanden | Experimental from a comparative study on correlation type registration algorithms[END_REF]. We briefly present the notations in Table 1 the five best measures of the different families that are considered.

(1) Family 1: Cross correlation-based measures -All these measures are based on a scalar product (Moravec, 1980) and NCC (Normalized Cross Correlation) is the most efficient one:

NCC(f l , f r ) = f l • f r f l f r . (1)
(2) Family 2: Classical statistics-based measures -These types of measures can be used: measures based on a distance or/and that are locally centered, variance-based or fourth-order cumulantbased measures [START_REF] Rziza | Dense disparity map estimation using cumulants[END_REF]. The best one is the LSAD (Locally scaled Sum of Absolute Differences) defined by:

LSAD(f l , f r ) = f l - f l f r f r 1 . (2)
(3) Family 3: Derivatives-based measures -Instead of using grey levels, these measures employ the derivatives of the images at different orders [START_REF] Seitz | Using local orientational information as image primitive for robust object recognition[END_REF]. Most of the existing measures use only the direction of the gradient vectors [START_REF] Ullah | Orientation Code Matching For Robust Object Search[END_REF], but, this kind of information can induce errors, in particular, with low norm gradient vectors whose direction is not reliable. In consequence, the most performant measure is based on the similarity of the image gradient vectors, GC (Gradient Correlation) [START_REF] Crouzil | A new correlation criterion based on gradient fields similarity[END_REF]. If the gradient vector at

p i, j w in I w is ∇I i, j
w and the norm is denoted by ∇I i, j w , the definition of GC is:

GC(f l , f r ) = ∑ A ∇I i+p, j+q l -∇I v+p,w+q r ∑ A ( ∇I i+p, j+q l + ∇I v+p,w+q r ) , (3) with
∑ A = ∑ N v p=-N v ∑ N h q=-N h .
(4) Family 4: Non-parametric statistics-based measures -They are based on the order of the grey levels inside the correlation window [START_REF] Kaneko | Robust image registration by increment sign correlation[END_REF][START_REF] Bhat | Ordinal measures for image correspondence[END_REF].

Using the order of the grey levels allows these measures to be robust against noises and occlusions but, sometimes, it also gives an ambiguous result, i.e. the best correlation score is obtained for the wrong correspondant. The most performant measure of this family is a non-parametric one, CENSUS [START_REF] Zabih | Non-parametric local transforms for computing visual correspondence[END_REF]. The similarity measure uses a transform that produces a bit chain which represents the pixels with an intensity lower than the central pixel:

R τ (f w ) = k∈[0;N f -1] ξ( f N f /2 w , f k w ),
where ξ( f

N f /2 w , f k w ) = 1 if f k w < f N f /2 w
and 0 elsewhere. CENSUS is the sum of the Hamming distances, denoted by D H , between the codes of each pixel of the correlation window:

CENSUS(f l , f r ) = N f -1 ∑ k=0 D H (R τ (f l ), R τ (f r )). (4)
(5) Family 5: Robust measures -We are particularly concerned with the occlusion problem which appears in the vicinity of a pixel near a depth discontinuity. In fact, some pixels lie on a first level of depth whereas the other pixels lie on a second level. It can disturb the matching process and introduce erroneous matches. To take this problem into account, robust statistics tools are introduced as correlation measures, like partial correlation [START_REF] Lan | Robust location based partial correlation[END_REF] or pseudonorms [START_REF] Delon | Analytic study of the stereoscopic correlation[END_REF]. The most efficient is SMAD, the Smooth Median Absolute Deviation [START_REF] Rousseeuw | L 1 -statistical analysis and related methods[END_REF]:

SMAD(f l , f r ) = h-1 ∑ k=0 (f l -f r -med(f l -f r )) 2 k:N f -1 ,
(5) where the ordered values f w are represented by:

( f w ) 0:N f -1 ≤ . . . ≤ ( f w ) N f -1:N f -1 .
It can be interpreted as a robust centered (median) and troncated distance and, in our experiments, h =

N f 2 .
Robust and non-parametric measures (families 4 and 5) are efficient in the presence of noises and/or occlusions whereas the classic ones (families 1 and 2) obtain better results when there is no major problems. The derivatives measures have been designed to be more efficient in the presence of noises, but, most of the time, they are really less efficient than the other ones, except GC which seems to have better results than the others, in particular in low textured areas. Interested readers can find more details about all the measures in [START_REF] Chambon | Occlusions handling in dense stereo matching[END_REF].

EVALUATION PROTOCOL

To validate our approach, 42 images, with their ground truth or reference disparity maps, have been tested (see Figure 1 for examples): 1 randomdot stereogram, 2 synthetic pairs (Murs) and one real image pair2 , and, finally, 38 real pairs introduced by Scharstein and Szeliski (9 in 2002 (Tsukuba) [START_REF] Scharstein | A taxomomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF], 2 in 2003 (Cones) [START_REF] Scharstein | High-Accuracy Stereo Depth Maps Using Structured Light[END_REF], 6 in 2005 and 22 in 2006 (Aloe)). The last ones are the most complex scenes. The most consequent evaluation protocol to highlight the different performances of global methods is given by the authors of [START_REF] Scharstein | A taxomomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF] 3 . Compared to their protocol, our comparison is based on all their 38 images instead of 4. about the estimation of these reference maps (ground truth), both in the cited papers and in the cited web page of section 3 (active vision is used and/or some constraints about the geometry of the scene are introduced).

Many criteria can be used to evaluate the quality of the results based on ground truth [START_REF] Chambon | Towards correlationbased matching algorithms that are robust near occlusions[END_REF].

However, this evaluation, we use the percentage of erroneous matches, noted ER, and the evaluation of the complementarity of the results (also based on ER) because they are the two most important aspects to consider in order to evaluate the impact of the proposed fusion algorithm.

COMPLEMENTARITY STUDY

To evaluate the complementarity of similarity measures, we analyse the percentage of erroneous matches (ER) for each measure used alone, and for each combination, by supposing that the correct correspondent is always kept (when one of the measures that are combined finds the exact correspondent), see Table 2 for the combination of 2 measures and Figure 2 for the percentage of erroneous matches with more than 2 measures. We use these notations: • M i , with i ∈ {1; . . . ; N m }, the N m tested measures;

• d th (p l ) the disparity of the pixel p l given by the ground truth; • d i (p l ), the disparity given by the algorithm based on the correlation measure M i ; • d N m tc (p l ) the theoretical or optimal combination of N m measures. More formally, the optimal combinations of the results over N m measures (N m ∈ {2; . . . ; 41} because in our previous work 41 measures have been studied), denoted d N m tc , is simply estimated by following this rule, for each pixel p l :

if

∃ i ∈ {1; . . . ; N m } where d i (p l ) = d th (p l ) then d N m tc = d th (p l ) and d N m tc is correct else d N m
tc is erroneous. We have tested these configurations: (C 1 ) N m = 2: All the 41 × 41 combinations have been evaluated and it highlights the best combination:

GC and SMAD with only 14.13% for the mean percentage ER on the 42 images.

(C 2 ) N m = 41: All the 41 measures have been theoretically combined and the results show that the percentage ER can be decreased to 7.26%.

(C 3 ) N m ∈ {3; . . . ; 40}: When we used the best combination GC-SMAD, any kind of measures can be added, the performances are quite equivalent. With 3 measures, the percentage ER decreases to about 13% and, then, it goes slowly to the minimum percentage ER (about 1% for each added measure) reached by the optimal combination of 41 measures. Moreover, when more than 10 measures are used, ER is close to this minimum.

First, the results show how the local matching with one correlation measure can be theoretically improved, and, second, which measures are the most complementary. In Table 2, we can remark that combining different measures highly improve the results: on the whole image, from 7% of improvement (2 measures combined) to 17% (41 measures). The following analyses illustrate these results. Using four maps, see Figure 3, we propose to visualize:

(1) The comparison of the two most complementary measures -As expected by the definitions of these measures, this visualization illustrates that SMAD compensates for the weaknesses of GC in occlusion areas (or near occlusion areas) whereas GC compensates for the weaknesses of SMAD in nonoccluded areas and, in particular, areas that are low textured, see (a) in Figure 3.

(2) The areas with 1 correct correspondent over 5, 10 or 41 correspondents -The most distinctive measure is GC, i.e. it is the most complementary measure to the other measures. SMAD is the second most complementary measure to the others, see (b) for 5, (c) for 10 and (d) for 41 in Figure 3. Moreover, with 10 different correlation measures, the results are quite near the results with 41 measures. And, our last conclusion is that combining more than 2 measures seems to be interesting because, most of the time, more than one measure obtains the correct correspondent. This last remark has inspired the fusion algorithm that is described in the next section. Figure 3: Study of complementarity (an example with image of Figure 1) -The pixels in grey levels correspond to pixels with more than 1 correct match over the N m combinations. The darker the pixel, the higher the number of correct matches. The color codes are used when only one measure gives the correct match. Moreover, we discern the occluded areas (Occ.) from the non-occluded areas (Non-occ.). In it shows that SMAD is efficient near occlusions whereas GC is more efficient than SMAD in low textured areas.

In conclusion, we have decided to present an algorithm that combines N different measures and we illustrate the interest of this kind of algorithm by using the two complementary measures: GC and SMAD.

In our first proposition of combination [START_REF] Chambon | Towards correlationbased matching algorithms that are robust near occlusions[END_REF], the algorithm was designed to take into account occlusions. In consequence, as expected, the results are good in non-occluded areas and also in occluded areas. The goal of this new algorithm is to improve this work by combining the advantages of each measure according to each kind of regions in order to take into account more difficulties, like low textured or noisy regions. Towards this goal, instead of detecting the occlusions, we work directly on how to merge the disparities by taking into account their variations in several matching maps (each map has been obtained with a different correlation measure).

Our method of fusion is based on two steps, the principle being to estimate a disparity map with each measure and then to merge the results applying the following two

(1) If more than one disparity map give the same match, the correspondence is validated and this result is considered as reliable.

(2) In an "undetermined area" (i.e. rule (1) is not respected), the "most reliable" disparity is kept.

The difficulty is to determine the most reliable. In this paper, we consider the disparities found in the neighborhood in the matching map of each considered measure.

Formally, these two rules can be defined as:

(1) Initialization for each pixel p l -The term d N m f is the final disparity, after the fusion of N m correlation measures.

If ∃ d | (d = argmax e M e (p l ) && (d ≥ N m 2 ) then d N m f (p l ) ← d else the disparity is undetermined. We define : M e (p l ) = #{i|d i (p l ) == e}.
(2) Refinement -For each pixel p l without disparity, we estimate the ambiguity, denoted by A, of each possible disparity d i (p l ). For the estimation of the function A, which represents how much the estimated disparity is reliable, we suppose that if most of the neighbors have the same disparity (in the same result obtained with the same correlation measure) the estimated disparity can be considered as sure. In consequence, for estimating A, we compare the studied disparity with the mean of the disparities in the neighborhood, denoted by N . The disparity with the lowest ambiguity is kept only if this ambiguity is not important, i.e. higher than a given threshold ε.

For (each pixel p l ) d = argmin i∈{1,N m } A(d i (p l )) with A(d i (p l )) = |d i (p l ) - 1 #N (p l ) ∑ k∈N (p l ) d i (p k )|
with N (p l ) the neighborhood of p l4 .

If (d < ε) 5 then d N m f (p l ) ← d else p l is occluded.

Matching results

For this part, the fusion algorithm has been tested with the fusion of the two most complementary measures: GC and SMAD. In order to try to detect occlusions and erroneous matches, we use the symmetry constraint that consists in estimating correspondences from the left image to the right image and then from the right to the left and in considering noncoherent matches as occluded pixels (these occluded pixels are shown in black in each disparity map). Table 3 shows the improvements of the percentage of erroneous matches obtained with the new algorithm of fusion on all the 42 tested images. The decreasing of this percentage is from 2.47 to 4.08 (with complex images), i.e. the images difficult to match because of the occlusion areas or the untextured areas. However, this improvement did not reach the theoretically maximal improvement that is showed in Table 2. Another way to appreciate the quality of the results is to look at the disparity maps that are given in Figure 4. The disparity maps obtained by fusion are the best ones because they contain less false negatives than the others. Moreover, the occlusion areas are better delimited (the contours are clean and contain no "holes").

As in the first step we have to estimate each disparity map induced by each measure, the execution time is the sum of the execution time of each correlationbased algorithm. The fusion algorithm does not take much time in comparison to the second step. In consequence, the higher the number of merged results, the higher the execution time and the execution time depends on the chosen measures. In our test, for example with Tsukuba, GC takes 17.6 s and SMAD 39.77 s, finally, the fusion algorithm takes about 1 minute. 4 for these images). The term Tc refers to the results obtained with a theoretical or optimal fusion, see Table 2. The percentage of erroneous matches with the new method is better than those obtained with the GC measure alone and in particular with complex scenes.

CONCLUSION

In this paper, we proposed a study of the complementarity of correlation measures, illustrated with visualization maps, and we introduced a new way to combine complementary measures. Moreover, we highlight the most complementary measures: GC and SMAD. The tests on 42 images illustrate the improvement of performances of the new fusion algorithm compared to classic correlation matching, i.e. based on one correlation measure alone. These results are encouraging but also exhibit the limit of this approach that might lead to investigate the fusion approach based on a voting method in the neighborhood of the studied pixel or to distinguish the most reliable measures (in the first step of the algorithm). Moreover, we will study the influence of the number of measures involved in the proposed algorithm.

  of data used in our tests (left images (a) and disparities 1 (b)). Interested readers can find more

Figure 2 :

 2 Figure 2: Percentage of erroneous matches versus the number of correlation measures theoretically combined. This graph illustrates the maximal number of measures that are interesting to combine (10) but it also highlights the biggest improvement obtained with only 2 measures.

Figure 4 :

 4 Figure 4: Disparity maps -(a), left image, (b) disparity map with SMAD, (c), with GC, (d), with FUSION. The fusion results present less false negatives, in particular for Cones and Aloe. The example of Tsukuba illustrates the limits of the method and the need to combine more than 2 measures.

Table 1 :

 1 Notations used for the description of the measures.

Table 2 :

 2

	MEASURE	ER	MEASURE	ER
	NCC	23.2	LSAD	23.3
	GC	21	CENSUS	20.2
	SMAD	27.9	GC-SMAD 14.13

Percentages of erroneous matches (ER) with each of the 5 best correlation measures and the best combinations of 2 correlation measures.

Table 3 :

 3 Precentage of erroneous matches -H represents the images with untextured areas, like Tsukuba pair, O, the images with a lot of occlusions, like Aloe pair and R, the images with no major difficulties, like Cones pair (seeFigure

A disparity (or matching) map represents for each pixel, the distance between the pixel and its correspondent. When the disparity map is represented by a grey level image, the clearer the pixel, the larger the distance is. Black pixels are occluded pixels.

http://www.irit.fr/˜Benoit.Bocquillon/MYCVR/ research.php

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

The 8 neighbors have been taken into account.

We have chosen ε = 1.