Vivien Maisonneuve
email: vivien.maisonneuve@cri.mines-paristech.fr.

Preservation of Lyapunov-Theoretic Proofs: From Real to Floating-Point Numbers

 Feron presents how Lyapunov-theoretic proofs of stability can be migrated toward computerreadable and verifiable certificates of control software behavior by relying of Floyd's and Hoare's proof system. However, Feron's proofs are addressed towards exact, real arithmetic and do not accurately represent the behavior of realistic programs run with machine arithmetic. We address the issue of errors resulting from the use of floating-point arithmetic: we present an approach to translate Feron's proof invariants on real arithmetic to similar invariants on floating-point numbers and show how our methodology applies to prove stability, thus allowing to verify whether the stability invariant still holds when the controller is implemented. We study in details the open-loop system of Feron's paper. We also use the same approach for Feron's closed-loop system, but the constraints are too tights to show stability in this second case: more leeway should be introduced in the proof on real numbers, otherwise the resulting system might be unstable.

Introduction

Provable stability constitutes an essential attribute of control systems, especially when human safety is involved, as in medical or aeronautical domains. Motivated by such applications, there exist many theorems to support system stability and performance under various assumptions and in various settings.

Stability criteria apply to a class of dynamical systems for which a stability proof is needed. Modern systems developments, such as adaptative control technologies, rely on robust stability and performance criteria as the primary justification for their relevance to safety-critical control applications. Lyapunov's stability theory plays a critical role in that regard.

The low-level software implementation of a control law can be inspected by analysis tools available to support the development of safety-critical computer programs. The simplest program analysis techniques consist of performing several simulations, sometimes including a software or hardware representation of the controlled system in the loop. However, simulations provide information about a large but only finite number of system behaviors. More advanced methods include model checking and abstract interpretation, e.g. using Astrée [START_REF] Cousot | The astrée static analyzer[END_REF]. In these methods, inputs are computer programs and outputs are certificates of proper program behavior along the chosen criterion. Another possibility is to use theorem-proving techniques, supported by computer tools such as Coq, Isabelle or PVS [START_REF]The coq proof assistant[END_REF][START_REF] Paulson | Isabelle[END_REF][START_REF] Owre | PVS specification and verification system[END_REF]. These proof assistants can be used to establish properties of programs and more general mathematical constructs. Model checking, abstract interpretation, and theorem-proving tools are all used to verify safety-critical applications.

In [START_REF] Feron | From control systems to control software[END_REF], Feron investigates how control-system domain knowledge and, in particular, Lyapunov-theoretic proofs of stability and performance, can be migrated toward computer-readable and verifiable certificates of control software behavior by relying on Floyd's and Hoare's proof system [START_REF] Peled | Software reliability methods[END_REF], applied to MATLAB pseudocode (Sections 2 and 3). His article focuses on using this framework to support such properties, namely, bounded-input, bounded-state stability, as they apply to control-system code implementations. But errors resulting from the use of floating-point arithmetic are not addressed. Feron presents both open-loop and closed-loop implementations of his controller, i.e. depending whether the system action on the environment and its feedback are modeled or not.

In this paper, we present an approach to translate Lyapunov-theoretic stability proof invariants on pseudocode with real arithmetic, as provided in Feron's article, to similar invariants on machine code that take into account rounding errors introduced by floating-point arithmetic. We use them to verify whether stability conditions still hold, in which case system stability with floating-point numbers can be established.

This document is organized as follows. The next section describes the second-order dynamical system example used by Feron in [START_REF] Feron | From control systems to control software[END_REF], with the corresponding controller. Next, Feron's analysis of the open-loop controller with real numbers is presented, followed by its rewriting with floating-point arithmetic. Then we discuss the case of the closed-loop system. The document concludes with a discussion of the generality of this successful approach.

We consider the first system described in the article of Feron [START_REF] Feron | From control systems to control software[END_REF]. It is a dynamical system composed of a single mass and a single spring shown in Figure 1. The position input y of the mass is available for feedback control. The signal y d is the reference signal, that is, the desired position to be followed by the mass.

A discrete-time MATLAB implementation of the controller, using real numbers, is provided in [START_REF] Feron | From control systems to control software[END_REF]. The source code is shown below: end Apart from the mechanical system state observation y and the desired system output y d , variables in this code are:

• x c = xc 1 xc 2 ∈ R 2
is the discrete-time controller state;

• y c ∈ [-1, 1] is the bounded output tracking error, i.e. the input y-y d passed through a saturation function to avoid variable overflow in the controller;

• u ∈ R is the mechanical system input, i.e. the action to be performed according to the controller.

Constants A c , B c , C c and D c are the discrete-time controller state, input, output and feedthrough matrices. The commands send and receive are basically I/O: they respectively send and receive data given in the commands first argument through a specific channel given by the commands second argument.

Open-Loop Stability Proof in Real by Feron

The stability proof of this system relies on Lyapunov theory. In simple terms, a system is Lyapunov stable if all states x c reachable from an initial starting state belonging to a bounded neighborhood V of an equilibrium point x e remain in V .

Lyapunov theory provides constraints that must be satisfied by such a V . On linear systems, they are equations that can be solved using linear matrix inequalities [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF]. Commonly, V is an ellipsoid.

In this case, to prove Lyapunov stability, we need to show that at any time, x c belongs to the set E P chosen by Feron according to Lyapunov's theory:

E P = {x ∈ R 2 | x T • P • x ≤ 1}, P = 10 -3   0.6742 0.0428 0.0428 2.4651   .
using E P as the stability neighborhood V . This set is a full ellipse, centered around 0 and slightly slanted:

x c ∈ E P ⇐⇒ 0.6742x 2 c 1 + 0.0856x c 1 x c 2 + 2.4651x 2 c 2 ≤ 1000.
E P is drawn in Figure 2.

A stability proof of the controller is provided in [START_REF] Feron | From control systems to control software[END_REF], using Floyd-Hoare program annotation technique [START_REF] Peled | Software reliability methods[END_REF]: each program instruction comes with an invariant. The program annotated by Feron is reproduced below.

xc = zeros(2,1); % x c ∈ E P receive(y, 2); receive(yd, 3); % x c ∈ E P while (1) % x c ∈ E P yc = max(min(y -yd, 1), -1); % x c ∈ E P , y 2 c ≤ 1 4 x c 1 -40 -30 -20 -10 10
; % ( xc yc ) ∈ E Qµ , Q µ = µP 0 2×1 0 1×2 1-µ , µ = 0.9991 u = Cc*xc + Dc*yc; % ( xc yc ) ∈ E Qµ , u 2 ≤ (C c D c ) • Q -1 µ • (C c D c ) -1 xc = Ac*xc + Bc*yc; % x c ∈ E P , P = (A c B c ) • Q -1 µ • (A c B c ) T -1 , u 2 ≤ (C c D c ) • Q -1 µ • (C c D c ) -1 send(u, 1); % x c ∈ E P receive(y, 2); % x c ∈ E P receive(yd, 3); % x c ∈ E P skip; % x c ∈ E P end
Most of the proof relies on algebraic arguments. For example, the invariant loosening on Line 9:

% x c ∈ E P , y 2 c ≤ 1 skip; % ( xc yc ) ∈ E Qµ , Q µ = µP 0 2×1 0 1×2 1-µ , µ = 0.9991 means x c ∈ E P ∧ y 2 c ≤ 1 =⇒ ( xc yc ) ∈ E Qµ with Q µ = µP 0 2×1 0 1×2 1-µ
and µ = 0.9991.

The correctness of this assertion stems from the fact that, given any value of µ ∈ [0, 1], the domain E Qµ is a solid ellipsoid, centered around 0, and whose intersection with the plane y c = 1 is equal to E P . Consequently, the solid bounded cylinder

C = {( xc yc ) | x c ∈ E P ∧ y 2 c ≤ 1} is included within E Qµ 5
(Figure 3).

E Qµ E P C x c 2 y c 1 -1
x c 1

Fig. 3. Inclusion of EP within EQ µ

The following invariants

% ( xc yc ) ∈ E Qµ , Q µ = µP 0 2×1 0 1×2 1-µ , µ = 0.9991 u = Cc*xc + Dc*yc; % ( xc yc ) ∈ E Qµ , u 2 ≤ (C c D c ) • Q -1 µ • (C c D c ) -1 xc = Ac*xc + Bc*yc; % x c ∈ E P , P = (A c B c ) • Q -1 µ • (A c B c ) T -1 , u 2 ≤ (C c D c ) • Q -1 µ • (C c D c
) -1 also rely on algebraic arguments and theorems.

Other invariants are trivial. Finally, only the very last loosening

% x c ∈ E P skip; % x c ∈ E P i.
e. E P ⊂ E P , that "closes" the loop, is not purely algebraic since its validity relies on the numerical parameters A c , B c , C c , D c . This assertion needs to be checked to ensure the correctness of the proof statements. This can be done either numerically, or algebraically for at most two-dimensional systems like this one. This stability proof has been checked with Mathematica [12] and is available online2 .

Open-Loop Stability Proof Scheme with Floating-Point Numbers

We would like to check that the stability proof still holds on a controller implemented with floating-point numbers. When using floating-point arithmetic, the values of constants are slightly altered and calculations are likely to produce rounding errors.

In absolute terms, it is impossible to switch from real to floating-point numbers without affecting the behavior of the controller. Thus, the stability proof cannot be preserved in the general case. On the other hand, the program still can be stable if rounding errors are small enough and the final inclusion E P ⊂ E P leaves them enough room. We study how proof invariants can be tweaked so that they correspond to a floating-point semantic. Our goal is to derive from the proof scheme for real numbers a proof scheme suited for floating-point arithmetic, whose correctness, although not guaranteed, can be checked as easily.

In this section, we use the IEEE 754 floating-point standard [6] encoded on 64 bits3 as most of today's floating-point units do. In this standard, both addition and multiplication are correctly rounded depending on the active rounding mode, which allows to bound the rounding error depending on the values of operands. The case of other numeric representations is discussed in Section 6.

Translation Scheme

The general translation scheme is the following (Figure 4). Each instruction instr R in the controller code with real numbers is turned into an equivalent instruction instr F on floating points, where real constants and operators are replaced by their floating-point counterparts. The post-condition post F on the floating-point side is computed using the same proof argument T as in the original code, enlarged to include the error term bound E that may occur in floating-point arithmetics in instr F . The precondition pre R might have been previously replaced by pre F as the postcondition of a preceding instruction.

The error term bound E depends on the instruction instr F and on the information on program variables given by pre F , hence the notation E(pre F , instr F ). The floating-point arithmetic error is not necessary bounded even when the variables are, e.g. in a division. In this case, we are not able to compute a suitable post F and the proof translation attempt fails.

Also, post F will be used as a precondition in the next instruction. An effort has to be made so that post F has a similar shape to post R , in order to fit with the proof argument T ′ used in this instruction, if possible. Thus, post F is an overapproximation of T(pre F , instr F ) + E(pre F , instr F ) in the domain of T ′ . Still, the analysis falls short if such a post F cannot be found.

% pre R instr R % post R = T(pre R , instr R ) % pre F instr F % post F ⊇ T(pre F , instr F ) + E(pre F , instr F )

Converting Constants

When translating Feron's program, the first step is to convert real constants: Apart from constants, the proof scheme for the first part of the program is unchanged:

xc = zeros(2,1); % x c ∈ E P receive(y, 2); receive(yd, 3); % x c ∈ E P while (1) % x c ∈ E P yc = max(min(y -yd, 1), -1); % x c ∈ E P , y 2 c ≤ 1 skip; % ( xc yc ) ∈ E Qµ , Q µ = µP 0 2×1 0 1×2 1-µ , µ = 0.9991 u = Cc*xc + Dc*yc;

Invariant on u

Then, we consider the next instruction in the original proof scheme:

% ( xc yc ) ∈ E Qµ , Q µ = µP 0 2×1 0 1×2 1-µ , µ = 0.9991 u = Cc*xc + Dc*yc; % ( xc yc ) ∈ E Qµ , u 2 ≤ (C c D c ) • Q -1 µ • (C c D c ) -1
First of all, matrices C c and D c must be replaced by their floating-point counterparts C f c and D f c both in the program instruction Line 10 and the ensuing invariant. This invariant relies only on algebraic arguments and does not depend on the values in the matrices, it still holds considering exact arithmetic operations. But this is not sufficient: indeed, this instruction is a sum of matrix multiplications, i.e. a set of additions and multiplications on floatingpoint numbers that yield rounding errors.

We can notice that entering this instruction, the values of matrices C f c , D f c and E Qµ are known, and the values of x c and y c are bounded by the precondition

( xc yc ) = xc yc ∈ E Qµ , that is 0.000673593x 2 c 1 + 0.000085523x c 1 x c 2 + 0.00246288x 2 c 2 + 0.9991y 2 c ≤ 1. (1)
From (1), we deduce:

           |x c 1 | ≤
|y c | ≤ 100 √ 9991 < 1.00046 (2)
Here we are able to find algebraic solutions, but this may be impossible with ellipsoids of higher dimension. Still, we would be able to find bounds using numerical methods.

In floating-point arithmetic, rounding errors created by addition and multiplication operators can be bounded when the operands are known or bounded by [START_REF] Cousot | The astrée static analyzer[END_REF], provided that overflow, underflow, and denormalized numbers do not occur [START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF][START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF].

Here, we need to compute

C f c x c + D f c y c = C f c (0,0) x c 1 + D f c y c 9
where all values in the right-hand term are known or bounded. Thus, a constant can be computed that bounds the absolute rounding error created when computing u. We used Rangelab [START_REF] Martel | Rangelab[END_REF], a static analysis tool to automatically validate the accuracy of floating-point or fixed-point computations, to compute an upper bound e for the error term4 : e = 5.90 • 10 -12 This way, starting from the algebraic result obtained on real numbers

|u| ≤ (C c D c ) • Q -1 µ • (C c D c ) T
we can ensure that with floating-point numbers, the inequality holds:

|u| ≤ U f = (C f c D f c ) • Q -1 µ • (C f c D f c ) T + e
which leads to the invariants:

% ( xc yc ) ∈ E Qµ , Q µ = µP 0 2×1 0 1×2 1-µ , µ = 0.9991 10 u = Ccf*xc + Dcf*yc; % ( xc yc ) ∈ E Qµ , u 2 ≤ U 2 f 4.4 Invariant on x c
The next instruction, considering constant changes, is:

11 xc = Acf*xc + Bcf*yc; % x c ∈ E Pf , Pf = (A f c B f c ) • Q -1 µ • (A f c B f c ) T -1 , u 2 ≤ (C c D c ) • Q -1 µ • (C c D c ) -1
where Pf is defined the same way P is, using floating-point terms A f c , B f c instead of the real-valued counterparts A c , B c . Again, this invariant holds independently of matrices values.

Here, we compute

   A f c (0,0) x c 1 + A f c (0,1) x c 2 + y c A f c (1,0) x c 1 + x c 2
Using the same method as above, absolute rounding errors introduced by floating-point operations can be bounded by constants These constants must be taken into account in the postcondition. Also, invariant on u must be replaced as previously. Then the postcondition can be replaced by

% x c ∈ E f Pf , u 2 ≤ U 2 f where E f
Pf is an ellipse that includes E Pf plus the rounding error terms (see Figure 5). Replacing the ellipse in the postcondition by another ellipse has the advantage of introducing little change in the stability proof sketch (instead of using a different domain, which would involve using different theorems), which can greatly facilitate tweaking the rest of the proof in longer codes. Formally, E f Pf must satisfy:

∀x c ∈ E Pf , ∀x ′ c ∈ R 2 , |x ′ c 1 -x c 1 | ≤ e 1 ∧ |x ′ c 2 -x c 2 | ≤ e 2 =⇒ x ′ c ∈ E f Pf (3)
E Pf E f Pf x c 1 x c 2 e 1 e 2

Fig. 5. Relation between E Pf and E f

Pf

At the end of the proof scheme, the system is stable with floating-point numbers if and only if the inclusion

E f Pf ⊂ E P holds. To succeed, E f
Pf should be as narrow as possible with respect to Equation (3). This is not a clear criterion, as several shapes are possible for E f Pf with no clear winner. We propose to define E f Pf as the smallest homothety of E Pf centered around 0 that satisfies (3). It can be computed rather easily, for any number of dimension; we give details for two dimensions.

Let a, b, c be the coefficients of E Pf :

E Pf = {(x c 1 , x c 2 ) | ax 2 c 1 + bx 2 c 2 + cx c 1 x c 2 ≤ 1}. 11
a, b and c are known, positive values. Then there exists k ≥ 0 s.t.

E f Pf = {(x c 1 , x c 2 ) | ax 2 c 1 + bx 2 c 2 + cx c 1 x c 2 ≤ k}. As E f Pf is wider than E Pf , k ≥ 1.
We need a condition on k that guarantees (3).

We consider a point (x c 1 , x c 2 ) located on the border of E Pf :

ax 2 c 1 + bx 2 c 2 + cx c 1 x c 2 = 1. (4)
By construction, for any values

ǫ 1 , ǫ 2 ∈ R s.t. |ǫ 1 | ≤ e 1 ∧ |ǫ 2 | ≤ e 2 , the relation
(x c 1 + ǫ 1 , x c 2 + ǫ 2 ) ∈ E f
Pf must hold, that is to say:

a(x c 1 + ǫ 1 ) 2 + b(x c 2 + ǫ 2 ) 2 + c(x c 1 + ǫ 1 )(x c 2 + ǫ 2 ) ≤ k.
It develops into

(ax 2 c 1 +bx 2 c 2 +cx c 1 x c 2 )+(2aǫ 1 +cǫ 2 )x c 1 +(2bǫ 2 +cǫ 1 )x c 2 +(aǫ 2 1 +bǫ 2 2 +cǫ 1 ǫ 2 ) ≤ k, that is 1 + (2aǫ 1 + cǫ 2 )x c 1 + (2bǫ 2 + cǫ 1 )x c 2 + (aǫ 2 1 + bǫ 2 2 + cǫ 1 ǫ 2 ) ≤ k.
due to [START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF].

Greatest values for the left-hand term are reached with ǫ 1 = |e 1 |∧ǫ 2 = |e 2 |, depending on the signs of x c 1 and x c 2 . As the ellipse E Pf is symmetric about the origin point (0, 0), we can set ǫ 1 = e 1 , which lets only two cases to study. Finally, we numerically verify that greatest values of the term are reached when ǫ 2 = e 2 . This is the only case we detail here.

We can write: 1 + αx c 1 + βx c 2 + γ ≤ k with values α = (2ae 1 + ce 2 ), β = (2be 2 + ce 1 ) and γ = (ae 2 1 + be 2 2 + ce 1 e 2 ). We known that x c 1 and x c 2 are bounded, thus so is the term αx c 1 + βx c 2 : we can compute a minimum bound δ s.t.

αx c 1 + βx c 2 ≤ δ. So it is sufficient that k satisfies: k ≥ 1 + γ + δ
Consequently, we define E f Pf as the smallest homothety of E Pf that satisfy (3), obtained with k = 1 + γ + δ. The instruction becomes:

% ( xc yc ) ∈ E Qµ , u 2 ≤ U 2 f xc = Ac*xc + Bc*yc; % x c ∈ E f Pf , u 2 ≤ U 2 f
In our case, starting from the ellipse that gives E f Pf .

E Pf = {(x c 1 , x c 2 ) | 0.00269007x 2 c 1 + 0.000341414x c 1 x c 2 + 0.

End of Proof Scheme

Then, what remains of the stability proof scheme becomes:

% x c ∈ E f P , u 2 ≤ U 2 f send(u, 1); % x c ∈ E f P receive(y, 2); % x c ∈ E f P receive(yd, 3); % x c ∈ E f P skip; % x c ∈ E P end
As previously, the final assertion E f Pf ⊂ E P must be checked. Two cases are possible:

• either E f Pf ⊂ E P , then we proved that the program is Lyapunov stable on a floating-point architecture;

• or E f In this open-loop case, we are able to check that E f Pf ⊂ E P (Mathematica script is available online: see footnote 2). Thus, the stability of the open-loop system with a 64-bit IEEE 754 compliant implementation is formally proven to hold using our proof translation scheme.

does not need to rely on the extensions of Hoare's logic to concurrent programs: one program at a time is running, through the blocking nature of the receive primitive.

Feron's stability proof with real numbers is much longer than for the openloop system. We do not detail it, the interested reader is referred to [START_REF] Feron | From control systems to control software[END_REF] for full information. To be noticed, the resulting comments are not much more complex than those available from the study of the controller alone. On the good side, as already mentioned, the Hoare formalism is not significantly affected by the concurrent structure of the closed-loop system.

A floating-point representation of the closed-loop system consists of keeping the right column of the listing above in its original settings, while replacing the left column with the corresponding floating-point implementation, as we did in Section 4. Using similar techniques to the study of the controller alone, proof invariants can be tweaked to take into account constant changes and rounding errors resulting from the use of floating-point arithmetic. Unfortunately, using these invariants it cannot be shown that the stability condition holds at the end of the loop body. In this case, we are unable to prove the system stability on a floating-point architecture: either the system is not stable with the floating-point based controller, or the proof parameters (E P , µ, . . . ) must be chosen more carefully by the controller designer.

Conclusion

The general idea is to replace some of the invariants in the original proof scheme by wider ones that include rounding errors, with the hope that the stability condition is strong enough and still holds. This approach is made possible by the fact that rounding errors introduced by the operations used in the code are bounded on bounded inputs and bounded controller state variables.

In this document, we study the case of a floating-point representation of real numbers. They are not available on all architectures, especially on microcontrollers that are commonly used to implement control systems. We quickly discuss the alternative real-number representations.

• We can deal with fixed-point arithmetic the same way we do with floatingpoint, as long as we stand far enough from extremal values that can lead to overflows.

• Another way to represent real numbers is to use two integers, a numerator and a denominator. Considering that the input values are exact, the elementary operations do not introduce rounding errors but can easily lead to overflows, e.g. when computing p 1 q 1 + p 2 q 2 = p 1 q 2 + p 2 q 1 q 1 q 2 .

A strategy must be used to prevent overflows by introducing approximations: in this case, the question is to quantify the errors introduced by these approximations.

In our example, we exclusively used additions and multiplications: divisions are not involved in linear control. Still, programs with divisions can also be analyzed, if the numerator can be shown to be far enough from zero: it is a supplementary constraint, but it is reasonable to assume that it should be respected on a realistic control system that uses divisions. Differentiable, periodic functions such as (sin) can be computed with an abacus and an interpolation function, thus with bounded error. In the same way, functions not periodic, but restricted to finite domains, can also be approximated. Other functions, such as (tan) or √ , should raise more issues.

Fig. 1 .

 1 Fig. 1. Mass-spring system

  (min(y -yd, 1), -1); skip; u = Cc*xc + Dc*yc; xc = Ac*xc + Bc*yc;

Fig. 2 .

 2 Fig. 2. The stability domain EP

Fig. 4 .

 4 Fig. 4. General translation scheme

  Ac = [0.4990, -0.0500; 0.0100, 1.0000]; Bc = [1; 0]; Cc = [564.48, 0]; Dc = -1280; to floating-point constants, assuming rounding to nearest value: Theses matrices A f c , B f c , C f c and D f c will be used instead of A c , B c , C c , D c in the sequel of the proof.

e 1 = 7 .

 17 42 • 10 -15 , e 2 = 3.62 • 10 -15

00247323x 2 c 2 ≤

 2 1}we get the following values:α = 1.03246 • 10 -17 , β = 1.84829 • 10 -17 , γ = 7.17582 • 10 -32 .Using Mathematica, we found δ = 5.35754 • 10 -16 ≫ γ and finally k = 1 + 5.35754• 10 -16

Pf⊂E

  P : as E f Pf was obtained through overapproximations, we cannot conclude about the program behavior.

Mathematica source file is available at: http://www.cri.ensmp.fr/people/ maisonneuve/lyafloat/resources/lyafloat_stability.nb, and the corresponding PDF file at: http://www.cri.ensmp.fr/people/maisonneuve/lyafloat/resources/ lyafloat_stability.pdf.

The procedure that follows would be exactly the same with 32-bit floating-point numbers, only with different numerical results.

Rangelab source file is available at: http://www.cri.ensmp.fr/people/maisonneuve/ lyafloat/resources/lyafloat_fp.m.

Closed-Loop Stability Proof

We now show how the proof of state boundedness of the closed-loop system specifications can be migrated to the level of the controller code and executable model of the system. To be more precise, we exploit the invariance of the ellipsoid E P to develop a proof of proper behavior, that is, stability and variable boundedness, for the computer program that implements the controller as it interacts with the physical system. Unlike the developments related to open-loop controller, this proof necessarily involves the presence of the physical system. In [START_REF] Feron | From control systems to control software[END_REF], Feron chooses to represent the physical system and the computer program by two concurrent programs, as shown below. In this scheme, the computer program representation of the physical system is to remain unchanged, since it only exists for modeling purposes and does not correspond to any actual program, whereas the controller code is allowed to evolve to reflect the various stages of its implementation.

Establishing proofs of stability of the closed-loop system at the code level is necessarily tied to understanding the joint behavior of the controller and the plant. The entire state space therefore consists of the direct sum of state spaces of the controller and the physical system. The approach described in the previous sections is used to document the corresponding system of two processes. One interesting aspect of these processes is their concurrency, which can complicate the structure of the state transitions. However, a close inspection of the programs reveals that the transition structure of the processes


Pauline Tan
email: pauline.tan@ens-cachan.fr

Pascal Monasse
email: monasse@imagine.enpc.fr

Stereo Disparity through Cost Aggregation with Guided Filter

Keywords: stereo-matching, cost volume, guided filter

Estimating the depth, or equivalently the disparity, of a stereo scene is a challenging problem in computer vision. The method proposed by Rhemann et al. in 2011 is based on a filtering of the cost volume, which gives for each pixel and for each hypothesized disparity a cost derived from pixel-by-pixel comparison. The filtering is performed by the guided filter proposed by He et al. in 2010. It computes a weighted local average of the costs. The weights are such that similar pixels tend to have similar costs. Eventually, a winner-take-all strategy selects the disparity with the minimal cost for each pixel. Non-consistent labels according to left-right consistency are rejected; a densification step can then be launched to fill the disparity map. The method can be used to solve other labeling problems (optical flow, segmentation) but this article focuses on the stereo matching problem.

Introduction

Stereo matching algorithms aim at estimating the depth of a scene given two photographs taken from different points of view (Figure 1). Depth estimation is done by estimating the disparity of each pixel, namely its apparent displacement from the first (reference) image to the second (target) image. This disparity permits to recover the 3D position of the point that was photographed in this pixel, if the stereo setup has been previously calibrated. We assume that the left image is the reference image. Let M be a 3D-point. We denote by O L (resp. O R ) the principal point of the left (resp. right) camera, that is the orthogonal projection of the optical center on the image plane. The point p (resp. q) denotes the image of M in the left image (resp. in the right image). Then, the disparity of M is given by d(M ) = --→ O R q ---→ O L p. Assume the camera makes a so-called fronto-parallel move, meaning that the focal plane is the same for all positions of the camera and that its motion is parallel to the x-axis of the camera CCD array. This is the case for instance in the Middlebury benchmark dataset 4 . Then, the disparity is scalar and its absolute value is inversely proportional to the distance from the observer. More precisely, if the camera moves from the left to the right, any point has a negative disparity (see Figure 2). Hence, the lower the (signed) disparity, the lower the depth.

In more general cases, the situation is slightly different. After stereo-rectification, some information is lost, such as the position of the camera principal points. The coordinates of pixels in the images are then not given with respect to the optical center, but to an arbitrary point (denoted O L and O R ). Without any supplementary information, the disparity of a point M is then given by d (M ) = --→ O R q ---→ O L p. We see that we can go back to the previous case by writing

d (M ) = ---→ O R O R + --→ O R q - --→ O L p - ---→ O L O L = d(M ) + C
where C does not depend on the point M . Then, to recover the depth from the disparity of any point, a constant has to be added to the computed disparity, but lower disparity still means lower depth.

If the camera motion is right to left, the situation is inverted and lower disparity means higher depth. This is important for the occlusion filling process, which in a first step fills in the non informed pixels with a disparity corresponding to lower depth.

L R × O L × O R • p • q • M Published in

Stereo Disparity through Cost A Filter

Pauline Tan 1 , Pascal M

1 CMLA, ENS Cachan, France (pauline 2 Imagine, ENPC, France (monasse@ L R

× O L × O R • p • q • M
Figure 2: Disparity when the camera moves from the left to the right. The disparity of the pixel M is given by subtracting the blue vector from the red vector. The disparity (in green) is then always negative.

Algorithm Description

Cost Volume

Pixel-to-pixel comparison is done by involving the truncated absolute difference of the color and the x-derivative of the grayscale images. These two terms are balanced with a parameter α.

1. Color comparison: for each pixel i in the left image and each pixel j in the right image, the color penalty is obtained by averaging the AD (Absolute Difference) on each color channel:

N color (i, j) := 1 3 I r L (i) -I r R (j) + I g L (i) -I g R (j) + I b L (i) -I b R (j) (2)
which is then thresholded:

N τ 1 color (i, j) := min N color (i, j), τ 1 . (3)
2.

x-derivative comparison: the discrete derivative along the x-axis of the gray-level image ĨL (resp. ĨR ), denoted ∇ x ĨL (resp. ∇ x ĨR ), is computed as follows5

∀ i = (i x , i y ), ∇ x ĨL (i) := ĨL (i x + 1, i y ) -ĨL (i x -1, i y ) 2 . ( 4
)
The AD is then computed

N gradient (i, j) := ∇ x ĨL (i) -∇ x ĨR (j) (5)
and thresholded N τ 2 gradient (i, j) := min N gradient (i, j), τ 2 . (6)

3. Eventually, the color and gradient terms are balanced:

N (i, j) := (1 -α) • N τ 1 color (i, j) + α • N τ 2 gradient (i, j). (7)
The cost volume C is then defined on a sampling integer grid of the volume [0,

M -1] × [0, N -1] × [d min , d max ],
C(i, d) := N (i, i + d). (8)
This assumes that i + d is in the domain of the right image (we note by a slight abuse of notation i

+ d := (i x + d, i y ) if i = (i x , i y )). If it is not the case, we set C(i, d) := (1 -α) • τ 1 + α • τ 2 which is the maximal possible cost.

Filtering with Guided Filter

The pixel-to-pixel comparison alone is too sensitive to noise. So, each slice p(i) = C(i, d) has to be filtered. The guided filter [START_REF] He | Guided image filtering[END_REF] is chosen since it is designed to preserve edges, by using the left image I L as guidance image.

Grayscale guidance image. Let I be the guidance (grayscale) image and q be the output of the filter. The key assumption of the guided filter is that it is a local linear model between I and q. For each square window ω k with radius r GF centered at k (by definition, a square with integer radius r is a square with side length 2 × r + 1 and |ω k | denotes the number of pixels in ω k ), a k-dependent filtered image q k is defined to be a linear transform of

I in ω k ∃ a k , b k , ∀i ∈ ω k , q k (i) = a k • I(i) + b k . ( 9
)
If k is too close to the image boundaries, then we only take into account the part of the window that remains in the domain of the image (see Figure 3). In the window ω k , since a k and b k are constant, we can write

∇q k = a k ∇I. (10)
That is, strong edges in I imply strong edges in q k . Hence, the guided filter preserves the edges of the guidance image. To determine coefficients a k and b k , the L 2 -norm between the filter input p and the filter output q k is minimized. Hence, an energy-based formulation of this problem is:

arg min a k ,b k E(a k , b k ) := i∈ω k a k • I(i) + b k -p(i) 2 + ε|ω k |a 2 k . ( 11
)
The first term is a data term and the second term is a regularization term preventing a k from being too large. The parameters a k and b k are then obtained by solving the Euler equation (∇E = 0):

a k = 1 |ω k | i∈ω k I(i)p(i) -µ k p k σ 2 k + ε (12) b k = p k -a k µ k ( 13
)
where µ k and σ 2 k are the mean and variance of I in the window ω k , and p k is the mean of p in ω k :

µ k := 1 |ω k | i∈ω k I(i), σ 2 k := 1 |ω k | i∈ω k I(i) 2 -µ 2 k and p k := 1 |ω k | i∈ω k p(i). (14)
The coefficients computed above depend on the window ω k . However, any pixel i is involved in several windows ω k , which leads to many values of coefficients a k and b k . A simple way to solve this problem is to take the average of all these values:

q(i) := 1 |ω i | k|i∈ω k q k (i) = 1 |ω i | k|i∈ω k a k • I(i) + b k = a i • I(i) + b i (15)
with a i (resp. b i ) the mean of a k (resp. b k ) in the window ω i , namely

a i := 1 |ω i | k|i∈ω k a k and b i := 1 |ω i | k|i∈ω k b k . ( 16
)
Notice that i ∈ ω k ⇔ k ∈ ω i , which explains the normalization factor 1/|ω i |.

The filter output q is then computed following four steps: Note that the first three steps use a box filter (mean in a square window). An efficient implementation of the box filter, independent of the radius of the square window, is shown in the next section. The above filter may be better understood when reformulated as a variable kernel convolution W GF :

q(i) = j W GF ij • p(j) with W GF ij := k∈ω i ∩ω j 1 |ω i | • |ω k | 1 + I(j) -µ k I(i) -µ k σ 2 k + ε . (17)
It is easily observed by this formula that pixels j in the neighborhood of i having similar color have a stronger weight in the resulting value of q(i).

Color guidance images. For color guidance images (I = (I r , I g , I b ) T ), the linear model can be written as

q k (i) = a T k I(i) + b k . (18)
The vectors a k and b k are then given by

a k = (Σ k + εI 3 ) -1          1 |ω k | i∈ω k I r (i)p(i) -µ r k p k 1 |ω k | i∈ω k I g (i)p(i) -µ g k p k 1 |ω k | i∈ω k I b (i)p(i) -µ b k p k          and b k = p k -a T k   µ r k µ g k µ b k   , (19)
with µ r k , µ g k , and µ b k the mean images of each color channel of the guidance image I:

µ [1] k := 1 |ω k | i∈ω k I [1] (i) for [1] ∈ {r, g, b}, (20)
Σ k is the covariance matrix of the guidance image I in the window ω k ,

Σ k :=   σ rr k σ rg k σ rb k σ rg k σ gg k σ gb k σ rb k σ gb k σ bb k   where σ [1][2] k := 1 |ω k | i∈ω k I [1] (i)I [2] (i) -µ [1] k µ [2] k for [1], [2] ∈ {r, g, b} (21)
and I 3 is the identity matrix of size 3 × 3. In this case, the kernel weights are given by

W GF ij = 1 |ω i | • |ω j | k∈ω i ∩ω j 1 + I(i) -µ k T Σ k + εI 3 -1 I(j) -µ k . (22)
In conclusion, each slice of the cost volume C(•, d) is filtered with the guidance image I L :

C filt (i, d) := j∈L 1 |ω i | • |ω j | k∈ω i ∩ω j 1 + I L (i) -µ k T Σ k + εI 3 -1 I L (j) -µ k • C(j, d). ( 23
)
The generalization to multi-channel images as guidance, as for example hyper-spectral satellite imagery, would follow the same process with matrices of different size.

Implementation of the Box Filter [3]

In this section, an implementation of the box filter proposed by Crow [START_REF] Crow | Summed-area tables for texture mapping[END_REF] is described. The simplest way to compute the mean in a square window of an image I is to sum the value of each pixel of the window, and then divide by the window size. Thus, for a window size of |ω k |, |ω k | -1 sums and one division are needed. Hence, the complexity of the process is linear with the radius of the window.

Another way to do it would be to first compute an integral image S of summed areas, of the same size as I:

∀ i = (i x , i y ), S(i) = ix kx=0 iy ky=0 I(k x , k y ). ( 24
)
Then a sum in any rectangle [i x , j x ] × [i y , j y ] is given by (see Figure 5)

S(j x , j y ) -S(i x -1, j y ) -S(j x , i y -1) + S(i x -1, i y -1), (25)
an eventual negative index removing the involved terms of S in the formula. Hence the computation of the mean only needs one sum, two subtractions and one division; it no longer depends on the radius of the window. 1 CMLA, ENS Cachan, France (pauline.tan@ens-cachan.fr) 2 Imagine, ENPC, France (monasse@imagine.enpc.fr)

0 0 i x j x i y j y

Disparity Selection

Once the cost volume is filtered, the disparity of any pixel i is the disparity that minimizes the cost associated to i in the cost volume ("Winner-take-all" strategy):

d(i) = arg min d ∈[d min ,dmax] C filt (i, d ) . ( 26
) (a)

Left-Right Consistency

The left-right consistency checks whether the disparity computed from the left image I L coincides with the disparity computed from the right image I R . At pixel precision, a disparity is rejected when

d i = -d i + d(i) (27)
where d denotes the disparity map computed from right to left image. This mask is expected to reject all occluded pixels and some flat region pixels of the scene. Thus rejected pixels are called "occluded pixels" in what follows. To handle subpixel disparities as well, this can be reformulated in the more general constraint:

d i + d i + d(i) ≥ d LR , (28)
with d LR = 1. Equation ( 27) is (28) when d LR = 0.

Densification Process

This step aims at filling the disparity map after the rejection test. Only occluded pixels are affected by this process, non-occluded ones keep their original disparity. It is done in two stages: the basic filling stage and the post-processing stage.

In the first stage a disparity value is given to the pixels rejected by the previous step. Since the mask mainly rejects occluded pixels, it often involves points of the scene that are hidden by other objects. In such cases, the rejected pixels are mainly in a region of discontinuity of the scene, where a closer object is in front of a farther one. Hence a strategy would be to adopt the disparity of the farther object, by exploring the closest nonoccluded pixels (on the same line). We recall that the farther object is given by the higher disparity. Hence, this leads to the following filling-in formula: Simple filling can generate streak-line artifacts. Hence an edge-preserving filter is used to remove them, by smoothing the disparity map. The left image is first filtered by the median filter (applied on 3 × 3-windows): for [START_REF] Boykov | A variable window approach to early vision[END_REF] ∈ {r, g, b},

∀ i = (i x , i y ), d fill (i) := max            d i x + min jx≤0 (ix+jx,iy ) nonoccluded {j x }, i y , d i x + min jx≥0 (ix+jx,iy ) nonoccluded {j x }, i y            . ( 29
) Notice that if i is nonoccluded, d fill (i) = d(i).
I [1] L,filt (i) = median I [1] L (i + j) j = (j x , j y ) : j x , j y ∈ {-1, 0, +1} . (30)
Let i be an occluded pixel. Then for each j in the square window ω i of radius r WMF centered on i, the bilateral filter weight is computed by

W BF ij := exp - i -j 2 σ 2 s • exp - I L,filt (i) -I L,filt (j) 2 c σ 2 c , (31)
where • denotes the Euclidean distance and • c the Euclidean distance in RGB color space:

I L,filt (i) -I L,filt (j) 2 c = I r L,filt (i) -I r L,filt (j) 2 + I g L,filt (i) -I g L,filt (j) 2 + I b L,filt (i) -I b L,filt (j) 2 . ( 32
)
This filter gives higher weight to pixels spatially close and of similar color, which are likely to belong to the same object. A disparity cumulative histogram h is then computed, where each value is weighted by the weight (31):

∀ d ∈ [d min , d max ], h(d) = j∈ω i |d(j)≤d W BF ij . (33)
The final disparity of i is then set to be the median value of h:

d final (i) := min d h(d) ≥ 1 2 h(d max ) . ( 34
)

Algorithm Summary and Implementation

A class Image defined in file image.h is used throughout. It is a container for an array of pixels that is reference counted: a copy of the image does not reallocate memory, it just adds a reference to the same array in memory (a coding pattern named shallow copy). The pixel by pixel operators of addition, subtraction and multiplication, used in the cost-volume filtering, are implemented as

Method Parameters

Several parameters should be tuned to use this method. Table 1 recalls all of them. Default values of these parameters suggested in [START_REF] Rhemann | Fast cost-volume filtering for visual correspondence and beyond[END_REF] are also given. Default values of τ 1 , τ 2 , ε and σ c are adapted to 8-bit images with intensity range in [0, 255] (most frequent case), but should be changed for higher bit depths.

Results obtained with different values of the parameters are then shown in order to see their influence. The experiments were done on the Tsukuba pair. Unless explicitly stated otherwise, only one parameter is changed at a time. The other parameters are tuned by default. Errors are also shown (absolute disparity error > 0.5) using a mask generated by Middlebury based on ground truth measurements (Figure 10). Note that this mask does not take into account the image borders. These results were also submitted to the Middlebury site in order to automatically get the percentage of bad pixels (with a wrong disparity). In the following figures (Figure 11 to Figure 18), the first row shows the disparity map obtained by our implementation. The second row shows the disparity map Compute the cost associated to pixel i for disparity d (7).

7

Set the cost value in p(i) = C(i, d) (8).

8

Cost-volume slice filtering by guided filter:

9
Compute the image pk of p (14).

10

Compute the images 1

|ω k | i∈ω k I r L (i)p(i), 1 |ω k | i∈ω k I b L (i)p(i), and
1 |ω k | i∈ω k I b L (i)p(i).
11

Compute the images of coefficients a k et b k (19).

12 Compute the images āi and bi (16).

13

Compute the output of the filter q = C filt (•, d) (15).

14

Disparity selection: Find the two nearest neighbors (on the same line) that are not occluded. Post-processing:

15 foreach pixel i in the left image do 16 if Cost(i)> q(i) then 17 Cost(i)← q(i) 18 d(i) ← d (26).
5

Filter the left image by median filter (30).

6

foreach occluded pixel i do 7

Compute the bilateral filter weights on a neighborhood of i (31).

8

Use these weights to compute a weighted median filter on the neighborhood of i (33), (34).

Parameter with the mask (in pink) generated by Middlebury. We first tested the influence of the parameters of the computation of the cost volume and its filtering by the guided filter. Figure 11 shows the results after the whole process (including filling and post-processing) when the size of the filter r GF is tuned. A larger filter size gives smoother results often visually more satisfying. However, the disparity map can be too smooth and errors generated on the boundaries of some objects. Indeed, when we look at the percentage of bad pixels on the boundaries of objects, we get 19.7% bad pixels for r GF = 19 while there are only 16.4% bad pixels for the default parameters. Nevertheless, the results remain globally better with a large filter size. Next, the influence of the parameter α was studied (see Figure 12). The larger this parameter, the heavier the gradient term weight. Results are better for large values of α. In other terms, the gradient comparison seems much more reliable than the color comparison. This is due to the fact that the scene is textured. Note that the disparity estimation is better on the lamp when α is small, because the lamp has no texture.

Figure 13 shows the error percentages for four of the Middlebury images with respect to r GF and α, the two main parameters and the ones that can be tuned in the online demo. This shows that the default parameters are well adapted to get optimal results in the Middlebury benchmark.

The other parameters needed for the computation of the cost volume are the thresholds τ 1 and τ 2 . When the filter size grows, the disparity map is cleaner (wrong disparities are removed, e.g. in the background), but errors appear on the object boundaries (e.g. on the space between the lamp and the table), due to an oversmooth map.

By default, they are chosen very low, which limits the influence of a large difference (of color or gradient). Figure 14 shows that the results are more accurate when these thresholds are larger. The results degrade when the thresholds are too large, though.

The last parameter involved in the filtering of the cost volume is the regularization parameter ε. The larger this parameter, the smoother the disparity map. Figure 15 shows that even though the number of bad pixels is lower for larger values of ε, the results are not visually satisfactory since the boundaries of many objects seem loose while ε grows.

The influence of the parameters involved in the filling process was also tested. We started with the size of the weighted median filter. Figure 16 shows that this parameter is not a key parameter, since the results do not significantly change when we tune it. The disparity estimation is slightly enhanced when the filter size grows. However, the larger the filter, the slower the computation.

Then, the parameters σ c (Figure 17) and σ s (Figure 18) are tuned, which control the importance of the color distance and the spatial distance in the bilateral filter. The experiments do not uncover any significant changes. One can notice that when these parameters are too large, the results are worse, because the bilateral filter is then not discriminatory enough. On the contrary, when they are too small, the filter rejects too many pixels.

By tuning properly the parameters to get better results, we managed to reach 9.01% bad pixels (and 8.44% in nonoccluded regions). Figure 19 shows the disparity map obtained with the following parameter values: r GF = 11, α = 0.95, τ 1 = 20, τ 2 = 10, ε = 255 2 • 10 -4 , r WMF = 29, σ c = 255 • 10 -1 , σ s = 15. However, one should keep in mind that the optimal values of the parameters depend on the scene.

We can compare the result before and after the densification process. Figure 20 shows the disparity map obtained with and without this step. The densification cleans the disparity map, by The larger the parameter, the better the results. Hence, the gradient comparison seems more reliable than the color comparison. This is specific to this scene (lots of textures). This is not true for scenes with many flat regions (the border of the lamp here is better estimated for a smaller value of α). Actually, by default, they are set too small. Thus, big differences do not affect the comparison too much. However, errors appear when they are too large (typically larger than those presented here). removing big errors (e.g. on the lamp). However, some details are lost (e.g. the wire on the lamp). Actually, one can notice that the densification step mainly enhances the result on flat regions of the disparity map. It does not affect much discontinuity regions. Then, one should consider this step to be a filtering step, as it removes big errors in flat regions. 5 Examples

(a) ε = 255 2 • 10 -4 (b) ε = 255 2 • 10 -1 (c) ε = 255 2 • 10 -2 (d) ε = 255 2 • 10 -5 (e) ε = 255 2 • 10 -4 (f) ε = 255 2 • 10 -1 (g) ε = 255 2 • 10 -2 (h) ε = 255 2 • 10 -5

Comparison with Original Algorithm

The authors of the original algorithm [START_REF] Rhemann | Fast cost-volume filtering for visual correspondence and beyond[END_REF] propose a Matlab implementation 6 . It performs similarly to our implementation, the slight differences being attributed to:

• Image color to gray conversion formulas.

• Handling of cost for pixels outside the image, term N (i, i + d) of ( 8): the Matlab code assumes a constant color R = G = B = b with b = 3 of I R outside its boundaries, while we assign the maximum possible cost.

• Numerical inaccuracies.

A comparison on the Middlebury stereo benchmark is given in Table 2. The authors' implementation is not their Matlab code. The latter uses d max = -1 instead of 0, but otherwise gives results

Error threshold=1: similar to our implementation. The CUDA code (not available) gives sensibly better results for the error threshold = 1 (only disparity errors strictly larger than 1 pixel are reported). Comparing both error thresholds, we see that our implementation has significantly more errors of 1 pixel. After a personal communication with the CUDA implementation author [START_REF] Hosni | [END_REF], it appears that on the Tsukuba pair, the left-right consistency threshold is d LR = 1 instead of 0. Whereas taking such a threshold yields slightly better results on Tsukuba and Cones, it deteriorates Venus and Teddy, so does not improve the final rank in the benchmark.

Moreover, we see in the green areas of Figure 21 some wrong disparity pixels that are not detected as occluded with the CUDA code but are still corrected in the final map. This is inconsistent with the described algorithm. No explanation could be found.

Results

We now show some other results obtained by our implementation. The experiments were done on images without ground truth. Hence, the results can only be appreciated visually. If not stated otherwise, the results shown were obtained by the whole process (including occlusion detection, filling and post-processing). The algorithm was launched with the default parameters (see previous section).

Globally, the results tend to show that this method gives very good results. The default parameters are fairly well chosen, which is sufficient most of the time. Hence, one can launch the code without tuning any parameter and yet get satisfactory results. Nevertheless, they can be enhanced by tuning the parameters, which are numerous. Thus, one has to tune and test each parameter several times to find the optimal parameters.

Chair. The code was tested on the Chair pair (Figure 22). Results are fairly good. However note that this pair is easy to deal with, since almost none of the stereovision difficulties are encountered (almost no reflection, textured regions). This result is compared with the result obtained before the left-right consistency check (Figure 23). Results are visually better. However, Figure 24 shows that this process degrades some parts of the disparity map, especially on some objects boundaries. Densification leads to better overall results, but some accuracy is lost.

Umbrella. Results on the Umbrella pair are shown in Figure 25. This experiment discloses a strong fattening effect, that is, neighbor points of objects in the foreground are labelled as if they belonged to the objects. Hence, such objects appear fatter than they actually are. This specific example shows that the boundaries of the green umbrella are not clearly drawn, since they are badly estimated (e.g. the left ear is bigger than expected). The background adopted the umbrella disparity.

Discussion on Adaptive-Window Approaches

The effectiveness of the method described in this article mainly relies on the use of the guided filter. Formula (17) expresses this filter with adaptive weights, which are shown in Figure 26. Hence, the filter only takes into account some selected neighboring pixels. The window is thus adapted to each pixel.

To compute the weights associated to a pixel i, we can use (17) or ( 23), but another possibility is for each pixel j to apply the guided filter to the Dirac image p j (k) = δ jk (1 at pixel j, 0 everywhere else) and take the filtered image at i, W GF ij = q j (i). Note that only pixels j at l ∞ -distance at most 2 • r GF from i should be considered, the sum in (23) reducing to 0 otherwise. Similar strategies, based on adaptive windows, have been proposed previously. Indeed, local stereo-matching methods highly rely on window comparison. Such comparison assumes that the depth within the window does not vary much, which is not the case on depth discontinuities. Hence, local methods require choosing the best window to perform well. The choice of a window depends on many criteria. Windows have to be large enough to avoid low signal to noise ratio. However, when they are too large, they may cover regions in which the depth varies a lot. Indeed, in such regions, the comparison may be affected by projective distortion due to strong depth discontinuities. Accordingly, several approaches have been tested to design adaptive windows. This strategy is aimed at selecting for each pixel the appropriate window.

One of the former approaches consists in adapting the window size to the pixel, knowing that it should depend on local variations of intensity and disparity [START_REF] Kanade | A stereo matching algorithm with an adaptive window: theory and experiment[END_REF]. Indeed, high variations of disparity mean depth discontinuities, which often entails distortions. Unfortunately, this method needs to evaluate the variation of the disparity, so it strongly depends on the initial disparity estimation.

Another approach is to adapt the shape of the window to fit with the shape of the discontinuities. Boykov et al. [START_REF] Boykov | A variable window approach to early vision[END_REF] find for each pixel an arbitrarily shaped connected window chosen by evaluating the plausibility of each pixel to have the same disparity as the pixel under consideration. Veksler [START_REF] Veksler | Stereo correspondence with compact windows via minimum ratio cycle[END_REF][START_REF]Fast variable window for stereo correspondence using integral images[END_REF] uses a set of windows (called compact windows) of different size and shape, and chooses the most reliable (with minimal average error). However, this method does not consider general shapes of window, which limits its performance.

More recently, a new adaptive-window method, proposed by [START_REF] Yoon | Adaptive support-weight approach for correspondence search[END_REF], uses square weighted windows. For each pixel, the weights are computed so that only neighboring pixels of interest (i.e. belonging to the same object) have heavy weight. Since the Gestalt principles state that grouping rules based on similarity and proximity are among the strongest, the weights should depend on the color similarity and the spatial proximity. The weights are chosen so that they decrease quickly while the pixels i

  (a) Left (reference) image (b) Right (target) image

Figure 1 :

 1 Figure 1: Tsukuba stereo pair

Let

  I L and I R be respectively the left and right color images of a stereo pair, of size M × N . Then, I r L , I g L and I b L denote the color channels of the left image (I r R , I g R and I b R for the right image). We proceed to the gray-level conversion of the images: ĨL = 0.299 • I r L + 0.587 • I g L + 0.0721 • I b L and ĨR = 0.299 • I r R + 0.587 • I g R + 0.0721 • I b R . (1) We assume the disparity range [d min , d max ] to be known.

k

  (a) Window of radius r GF = 2 centered at k = (3, 6). |ω k | = 25.

k

  (b) Window of radius r GF = 2 centered at k = (1, 8). |ω k | = 16.

Figure 3 :

 3 Figure 3: Examples of square windows.

  (a) Guidance image (b) Slice d = -4 of the cost volume (c) Filtering by guided filter

Figure 4 :

 4 Figure 4: Filtering by guided filter of a slice of the cost volume for the Tsukuba pair (τ 1 = 7, τ 2 = 2, α = 0.9, ε = 10 -4 , r GF = 9).

Figure 5 :

 5 Figure 5: Calculation of a summed area from the integral image.

  Left image (b) Right image (c) Disparity map

Figure 6 :

 6 Figure 6: Result of the winner-take-all algorithm on the Tsukuba pair.

Figure 7 :

 7 Figure 7: Result of the left-right consistency check on the Tsukuba pair. Most of the occluded pixels are detected by the left-right consistency mask. So are some pixels on constant disparity regions.

  (a) Left image (b) Left-right consistency rejections (c) Simple filling

Figure 8 :

 8 Figure 8: Result of the filling on the Tsukuba pair. Streak-line artifacts appear.

Figure 9 :

 9 Figure 9: Result of the post-processing on the Tsukuba pair. The disparity map is smoothed.

Algorithm 1 : 1 Cost(.)← ∞ 2 foreach d ∈ [d min , d max ] do 3 Building the slice d of cost volume: 4 foreach pixel i in the left image do 5

 112345 Cost-volume guided filtering algorithm Input: Color images I L and I R , disparity range [d min , d max ]. Output: Disparity map d(.) Compute the images µ r k , µ g k , µ b k (20), and Σ k (21).Compute the color term (3) and the gradient term (6) for disparity d.

6

 6

19 Algorithm 2 : 1 foreach pixel i in the left image do 2 3 Algorithm 3 : 1 foreach occluded pixel i do 2

 2123312 Left-right consistency check Compute the disparity map d from I R to I L . If |d(i) + d (i + d(i))| > d LR , label the pixel i as occluded. Disparity map filling Simple filling:

3

  Adopt the higher (assuming left to right camera motion) disparity (29).

4

 4

  (a) Left image (b) Ground truth

Figure 10 :

 10 Figure 10: Ground truth used to evaluate the results.

Figure 11 :

 11 Figure 11: Influence of the guided filter size. (a)+(e) Default parameters (r GF = 9). 12.9% are bad pixels. (b)+(f) r GF = 4, bad pixels: 15.9%. (c)+(g) r GF = 11, bad pixels: 11.8%. (d)+(h) r GF = 19, bad pixels: 11.7%.When the filter size grows, the disparity map is cleaner (wrong disparities are removed, e.g. in the background), but errors appear on the object boundaries (e.g. on the space between the lamp and the table), due to an oversmooth map.

1 Figure 12 :

 112 Figure 12: Influence of the parameter α. (a)+(e) Default parameters (α = 0.9), bad pixels: 12.9%. (b)+(f) α = 0.95, bad pixels: 10.7%. (c)+(g) α = 0.5, bad pixels: 20.9%. (d)+(h) α = 0.1, bad pixels: 23.5%.The larger the parameter, the better the results. Hence, the gradient comparison seems more reliable than the color comparison. This is specific to this scene (lots of textures). This is not true for scenes with many flat regions (the border of the lamp here is better estimated for a smaller value of α).

Figure 13 :

 13 Figure 13: Influence of the parameters r GF (left) and α (right) on errors in the Middlebury benchmark.

(a) τ 1

 1 = 7, τ 2 = 2 (b) τ 1 = 20, τ 2 = 2 (c) τ 1 = 7, τ 2 = 10 (d) τ 1 = 20, τ 2 = 10 (e) τ 1 = 7, τ 2 = 2 (f) τ 1 = 20, τ 2 = 2 (g) τ 1 = 7, τ 2 = 10 (h) τ 1 = 20, τ 2 = 10

Figure 14 :

 14 Figure 14: Influence of thresholds τ 1 and τ 2 . (a)+(e) Default parameters (τ 1 = 7, τ 2 = 2), bad pixels: 12.9%. (b)+(f) τ 1 = 20, bad pixels: 12.1%. (c)+(g) τ 2 = 10, bad pixels: 9.65%. (d)+(h) τ 1 = 20, τ 2 = 10, bad pixels: 9.49%. The results improve when the thresholds are chosen a little larger.Actually, by default, they are set too small. Thus, big differences do not affect the comparison too much. However, errors appear when they are too large (typically larger than those presented here).

Figure 15 :

 15 Figure 15: Influence of the regularization parameter. (a)+(e) Default parameters (ε = 255 2 • 10 -4 ). Bad pixels: 12.9%. (b)+(f) ε = 255 2 • 10 -1 , bad pixels: 13.5%. (c)+(g) ε = 255 2 • 10 -2 , bad pixels: 12.2%. (d)+(h) ε = 255 2 • 10 -5 , bad pixels: 13.0%. For large values of ε, the object boundaries are loose and the disparity map is smoother.

  Figure 16: Influence of the size of the weighted median filter. (a)+(e) Default parameters (r WMF = 9), bad pixels: 13.0%. (b)+(f) r WMF = 4, bad pixels: 13.2%. (c)+(g) r WMF = 19, bad pixels: 12.9%. (d)+(h) r WMF = 29, bad pixels: 12.8%. This parameter shows no influence on the results. However, as expected, the results are slightly better when the filter size grows.

4 Figure 17 :

 417 Figure 17: Influence of the parameter σ c . (a)+(e) Default parameters (σ c = 255 • 10 -1 ), bad pixels: 12.9%. (b)+(f) σ c = 255 • 10 -2 , bad pixels: 13.3%. (c)+(g) σ c = 255, bad pixels: 12.7%. (d)+(h) σ c = 255 × 4, bad pixels: 12.8%. This parameter has no strong influence, but it should not be chosen too large or too small.

  Figure 18: Influence of the parameter σ s . (a)+(e) Default parameters (σ s = 9), bad pixels: 12.9%. (b)+(f) σ s = 10 -1 , bad pixels: 13.4%. (c)+(g) σ s = 2, bad pixels: 13.3%. (d)+(h) σ s = 15, bad pixels: 12.7%. This parameter does not change the results a lot, but it should not be chosen too large or too small.

Figure 19 :

 19 Figure 19: Enhancement by tuning properly the parameters. (a)+(c) Default parameters. Bad pixels: 12.9%. (b)+(d) Tuned parameters.Bad pixels: 9.01%. We managed to enhance the result, especially in the background of the scene, while keeping a good estimation anywhere else. Nevertheless, one can notice that the disparity map is smoother, which creates some errors (e.g. on the right border of the lamp or on the statue).

Figure 20 :

 20 Figure 20: Enhancement by the densification process. (a)+(c) Without densification. Bad pixels: 14.5%. (b)+(d) With densification. Bad pixels: 12.9%. The densification mainly enhanced the results on nonoccluded regions. It removed some obvious errors easily detected (e.g. on the lamp or on the statue), but it also removed some details (e.g. the wire on the lamp).

Figure 21 :

 21 Figure 21: CUDA implementation result on Tsukuba pair. Left: disparity map with occluded pixels. Right: final disparity map.

Figure 22 :

 22 Figure 22: Chair pair. (a) Left image. (b) Right image. (c) Result. The result is very satisfactory, the depth of the scene seems probable.

Figure 23 :

 23 Figure 23: Chair pair. (a) Left image. (b) Before the occlusion detection. (c) After the whole process. After the occlusion detection and the densification step, the disparity map is cleaner, most of obvious errors have been removed. Moreover, the disparity has been well interpolated on the left border of the image.

Figure 24 :

 24 Figure 24: Chair pair (detail). (a) Left image. (b) Before the occlusion detection. (c) After the whole process.In this detail, the disparity was better estimated before the occlusion detection, even though the densification process removed some errors.

Figure 25 :

 25 Figure 25: Umbrella pair. (a) Left image. (b) Right image. (c) Result.

Figure 26 :

 26 Figure 26: Adaptive weights in the guided filter (left image for each patch) and of Yoon-Kweon's bilateral filter [14] (right image for each patch).

  such that, for any disparity value d ∈ [d min , d max ]:

Table 1 :

 1 Method parameters.

		Suggested value Description
	r GF	9	Guided filter radius.
	α	0.9	Parameter in [0, 1] that balances the color and the
			gradient contribution. The bigger this parameter,
			the more important the contribution of the gradi-
			ent term.
	τ 1	7	Threshold used for the color penalty AD distance.
	τ 2	2	Threshold used for the gradient penalty AD dis-
			tance.
	ε d LR	255 2 • 10 -4 0	Guided filter regularization parameter. Left-right disparity difference tolerance
	r WMF	9	Weighted median filter radius.
	σ s	9	Weight of the spatial term in the bilateral filter.
	σ c	255 • 10 -1	Weight of the color term in the bilateral filter.

Table 2 :

 2 Results on Middlebury stereo benchmark (AR=Average Rank). First implementation is our code, the second one is the authors' CUDA code results, as stored in the benchmark.

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

In[START_REF] He | Guided image filtering[END_REF], pixel intensities outside the image are considered to be the same as at the closest pixel inside the image.

https://www.ims.tuwien.ac.at/publications/tuw-202088

The original article[START_REF] Yoon | Adaptive support-weight approach for correspondence search[END_REF] uses Euclidean distance in CIELab color space, in our experiments we use L 1 distance in RGB space.

and j are dissimilar or not close enough:

where • c denotes the Euclidean distance between two colors represented in color space 7 . Note that this is close to the bilateral filter used in the post-processing. Figure 26 shows examples of weights obtained by this formula. A limitation of this method is that it is computationally expensive, while the guided filter discloses similar results with a lower complexity. Note that the name "bilateral filter" [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF] is the most recent one for a filter known as the sigma filter [START_REF] Lee | Digital image smoothing and the sigma filter[END_REF], the Yaroslavsky filter [START_REF] Yaroslavsky | Digital picture processing: an introduction[END_REF] and the SUSAN filter [START_REF] Smith | Susan-a new approach to low level image processing[END_REF]. The bilateral weights in Figure 26 are computed with parameters: σ 2 s = 17.5 (half the window size 35 × 35) and σ 2 c = 14 • 3 (3 channels). The weights appear roughly equivalent by both methods. Note that there is implicitly a decreasing tendency with the distance to central pixel i in the guided filter, because the further pixel j, the fewer windows ω k containing i and j. A full comparison with other adaptive-window methods will be the object of future work.

Image Credits

All images by the authors (license CC-BY-SA) except:

Middlebury.