Marcin Konczykowski

Cornelis Jacominus Van Der Beek

Alexei E Koshelev

Vincent Mosser

Matthew Dodgson

Pieter H Kes

C J Van Der Beek

Composite to tilted vortex lattice transition in Bi_2Sr_2CaCu_2O_8 + δ in oblique fields

Keywords: numbers: 74.25.Qt, 74.25.Op, 74.25.Dw

come

The first order "vortex melting" transition from a solid (phase-ordered) state to a liquid state with only short range correlations is the main feature of the phase diagram of vortex lines in clean, layered high-temperature superconductors [1]. The application of a small field component H , parallel to the superconducting layers, leads to a lattice of tilted vortex lines that melts in a similar fashion [2]. However, in the more anisotropic (layered) compounds such as Bi 2 Sr 2 CaCu 2 O 8+δ , the depression of the perpendicular component of the melting field H ⊥ m by larger parallel fields was interpreted as the consequence of the decomposition of the tilted vortex lattice into a combined lattice structure of Josephson Vortices (JVs) and Abrikosov-type pancake vortices (PVs) [2]. For very small field components H ⊥ perpendicular to the layers, chain structures [3] arising from the attractive interaction of PVs with JVs were directly visualized by Bitter decoration [4,5], scanning Hall-probe [3,6] and magnetooptical techniques [7,8]. At higher H ⊥ ∼ H ⊥ m , the contribution of the JVs to the free energy of the pancake vortex crystal results in the almost linear depression of H ⊥ m as function of the parallel field [2,10,11,12]. This behavior in moderate H stops at a temperature dependent characteristic field H cr . Even though melting is still observed above H cr , the variation of H ⊥ m with increasing H becomes much weaker [11,12]. Several controversial interpretations of this changing behavior were proposed, such as layer decoupling [11], a commensurate transition [13], and a matching effect [14].

In this Letter we focus on the high-temperature portion of the vortex phase diagram in single crystalline Bi 2 Sr 2 CaCu 2 O 8+δ in oblique fields, which can be established precisely using the well-defined discontinuity of the vortex density at the melting transition. We show that (H ⊥ m , H cr ) corresponds to a tricritical point in the vortex lattice phase diagram, where the melting crosses a novel transition from a composite lattice at low parallel fields, to another tilted lattice structure at high H . The experimental observation of large hysteresis suggests that this transition is strongly first order, consistent with recent predictions [15]. The identification of the vortex ground state at high parallel field as a tilted lattice structure resolves the open problem of the apparent anisotropy factor γ ef f , and allows one to determine the enhancement of H ⊥ m by magnetic coupling. We find the temperature dependence of γ ef f to be consistent with previous observations [16,17] and in quantitative agreement with the proposed model.

Experiments were performed on rectangular samples cut from Bi 2 Sr 2 CaCu 2 O 8+δ single crystals with different oxygen content [18]. The c-axis component of the local magnetic induction B ⊥ (r) was measured by micro-Hall sensors placed on the central part of the sample. The 2D electron gas Hall sensors were fabricated in GaAlAs heterostructures and had 8 × 8 µm 2 active area. Results are presented in Fig. 1(a) as the local magnetization H ⊥ s ≡ B ⊥ -H ⊥ . The local dc magnetization of all crystals shows a sharp discontinuity, ∆B ⊥ , at the vortex melting transition, that was tracked as function of H at various fixed temperatures. The angle θ between the magnetic field and the crystalline c-axis was computercontrolled with 0.001 • resolution, while the field magnitude could be swept up to 1 T using an electromagnet. Two types of magnetization loops were measured. In the first, the magnetization is traced as function of the c-axis field at constant H ; in the second, the magnetization is measured as function of H at constant H ⊥ . While the discontinuity in the dc magnetization gives a clear identification of the melting field, another method [19], in which the magnitude B(f, T ) of the periodic part of the induction above the sample is measured at the frequency f of an ac ripple field applied perpendicularly to the sample plane, is more convenient and precise. The ac response is represented as the transmittivity T ′ , i.e.

T ′ ≡ [B ′ (f, T ) -B(f, T ≪ Tc)]/[B(f, T > Tc) -B(f, T ≪ Tc)].
the in-phase component B ′ (f, T ), normalized by the amplitude h ac of the ac ripple [21]. The steplike feature in the dc magnetization loop at H ⊥ m translates to a paramagnetic peak in the ac response, shown in Fig. 1(b) [19]. The magnitude of this peak depends on the ratio of ∆B ⊥ to h ac . The peak position is independent of both the amplitude and frequency of the ac ripple. In the explored temperature range (above 50 K) and at low frequency (below 27 Hz), a true paramagnetic signal is measured. At higher frequencies or lower temperatures, flux pinning results in the partial shielding of the ac field [20]. Nevertheless, a peak-like feature persists at melting.

Figure 1(a) shows that at T > 50 K, the application of even a small magnetic field component parallel to the layers results in the drastic suppression of magnetic irre- versibility. This is expected when the geometric barrier is at the origin of flux pinning [22,23]. Simultaneously, H ⊥ m is depressed linearly with increasing H . However, at a well-defined value H cr , the dependence of H ⊥ m on in-plane field changes to a much slower, quadratic behavior that very well fits the anisotropic London model, [24]; i.e. the perpendicular component of the melting field

H m (θ) = H m0 /(cos 2 θ + sin 2 θ/γ 2 ef f ) 1/2
H ⊥ m = H ⊥2 m0 - H 2 γ 2 ef f ≈ H ⊥ m0 1 - H 2 2γ 2 ef f H ⊥2 m0 . (1)
The characteristic field H m0 and the effective anisotropy parameter γ ef f will be defined below. In Fig. 1(b), another feature in the in-phase component of the ac response can be distinguished, at perpendicular fields H ⊥ somewhat smaller than H ⊥ m . This feature is brought out much more clearly in sweeps of the parallel field, shown in Fig. 2. There is an abrupt jump from lower to higher values of T ′ on increasing H , that only appears for parallel fields H < ∼ H cr . The position of the jump does not depend on ac frequency. At low amplitude of the ac field, a pronounced hysteresis of T ′ is observed; this disappears if h ac is sufficiently increased.

The transmittivity T ′ is simply related to the magnitude of the shielding current flowing in the sample in response to the applied ac magnetic field, a higher T ′ corresponding to a smaller current and less screening [21]. In the present case, dc magnetization loops point to the geometrical barrier [22] as the main source of screening. However, increasing the ac field frequency reduces the role of thermally activated depinning of vortices in the crystal bulk; as a consequence, a bulk screening current with the first-order melting transition H ⊥ m (•) determined from the paramagnetic peak in T ′ , and the first-order transition to the tilted PV lattice at Bt (•), determined from the "glitch" in T ′ (Fig. 1 b). The dashed line is a fit to the composite-to-tilted lattice transition, Eq. ( 2) with C = 0.030; the continuous line is a fit of the high-field portion of the vortex lattice melting line to Eq. ( 1).

due to vortex pinning emerges [20]. At the frequencies of Fig. 2, the step in T ′ is due to a discontinuous change of the magnitude of this bulk current at the well-defined in-plane field, H ≡ H ct . The location of H ct does not depend on the frequency and h ac which indicates a vortex phase transition in the bulk, from a low H -phase with higher pinning, to a high H -phase with lower pinning. The hysteresis of the screening current indicates it to be first order.

In Fig. 3 we collect, for T = 75 K, the positions of the two first order transitions in a plot of H ⊥ versus H . The usual melting field H ⊥ m of the vortices, deduced from the paramagnetic peak in the transmittivity, shows the well-known linear decrease as function of H [2,10,11,12], up to the field component H cr . The field H ⊥ ct of the first order transition revealed by the (irreversible) transmittivity rapidly increases with H and crosses the melting line at H cr in a tricritical point. The same scenario is observed at all explored temperatures (T > 50 K), with temperature dependent values of H cr . The anisotropy factor γ ef f , extracted from the London model fit to the high-H part of the melting line, depends on temperature as well as on the oxygen content of the Bi 2 Sr 2 CaCu 2 O 8+δ crystals; it is depicted in Fig 4.

The low-H portion of the phase diagram with almost linear H ⊥ m (H ) dependence has been interpreted as the region of crossing vortex lattices of JVs and PV stacks [2]. A more accurate analysis shows that other "compositelattice" configurations compete for the ground state in the parameter range of Bi 2 Sr 2 CaCu 2 O 8+δ . These are the soliton lattice [25], as well as the set of combined lattices composed of regularly spaced rows of tilted pancake stacks, separated by M rows of pancake stacks aligned with the c-axis. The latter type of lattice becomes favorable at smaller anisotropies and larger H . Moreover, if H is sufficiently large and the material anisotropy is not extremely high, a simple tilted lattice (M = 0) turns out to be the most favorable configuration. We interpret the experimentally observed transition as that from a composite to such a uniformly tilted lattice. A simple estimate for the in-plane field at which this transition is expected, B ct , can be obtained by comparing the ground state energies of the simplest (M = 1) composite lattice and of the uniformly tilted lattice, giving [15]

B ct ≈ C γ λ ab B ⊥ Φ 0 / ln 1.55 B ⊥ Φ 0 sB ct 1/2 . ( 2
)
Here Φ 0 is the flux quantum, λ ab is the ab-plane penetration depth, γ is the penetration depth ratio λ c /λ ab , and s is the layer spacing. Eq. ( 2) gives a very good fit to the experimental transition line, as illustrated in Fig. 3. The anisotropic three-dimensional behavior (1) of H ⊥ m for large in-plane field B > B cr strongly supports this interpretation. The H ⊥ m (H )-dependence is the direct consequence of the vanishing contribution of the magnetic interaction between PVs to the vortex tilt stiffness in a highly inclined tilted vortex structure. The angular dependence of the melting field can be derived using a scaling transformation of coordinates, z = γ 2/3 z; r⊥ = γ -1/3 r ⊥ , which reduces the larger part of the free energy to an isotropic form [24]. In scaled coordinates the magnetic field is given by B = Bγ 2/3 cos 2 θ + sin 2 θ/γ 2 1/2 , while the tilt angle tan θ = tan θ/γ. The Josephson tilt energy of a deformed vortex line (PV stack) in scaled coordinates,

E J,t = d kl 2π ε1 ( kl ) 2 k2 l |δũ( kl )| 2 ,
is determined by the effective line tension ε1 ( kl ) = ε0 ln 1/r cut kl , valid when the wave vector along the line direction, kl , is much larger than the vortex lattice zone boundary vector. Here, δũ( kl ) is the Fourier transform of the line deformation, and ε0 = ε 0 γ -2/3 with ε 0 ≡ Φ 2 0 /(4πλ ab ) 2 . For near-perpendicular fields ( θ ≪ 1) the core cut-off distance rcut is determined by the so-called thermal vortex wandering length, rcut ≈ ũ2 n,n+1

1/2 ≡ (ũ n+1 -ũn ) 2 1/2
, where u n is the position of the PV vortex in layer n [26]. For a tilted vortex line, ũn,n+1 = s tan θ + δũ n,n+1 consists of the average displacement as well as random (thermal) fluctuations meaning that the core cut-off r2 cut ≈ s 2 tan 2 θ + (δũ n,n+1 ) 2 . The melting temperature is given by T m = A ε1 (1/ã)ε 0 ã, with ã ≡ (Φ 0 / B) 1/2 and A ≈ 0.1 [27]. Returning to real coor- FIG. 4: Temperature dependence of the apparent anisotropy γ ef f extracted from the fits of the melting line to Eq. ( 1). The drawn line shows a fit to Eq. ( 4) with intrinsic γ = 500.

dinates, we obtain

T 2 m = A 2 (ε 0 s) 2 ln C J B sc (θ)/B r 2 0 + tan 2 θ/γ 2 B sc (θ) B (3)
where the numerical constant C J ≈ 5 can be estimated within the self-consistent harmonic approximation, B sc (θ) = (Φ 0 /γ 2 s 2 )/ cos 2 θ + γ -2 sin 2 θ, and

r 2 0 = (δũ n,n+1 ) 2 /(γs) 2 ≈ 2A[Φ 0 /s 2 γ 2 B m (0)] 1/2
. Note that the angular-dependent core cutoff introduces an additional angular dependence of melting field: T m no longer depends only on the ratio B sc (θ)/B. In particular, a new angular scale appears given by tan θ = γr 0 . In the experimental angular range tan θ ≪ γ, γr 0 , we recover Eq. ( 1) with the apparent anisotropy

γ ef f ≈ γ   1 + 10 B m (0)γ 2 s 2 /Φ 0 ln 68 Φ 0 /(B m (0)γ 2 s 2 )   -1/2 . (4)
We note several key points. First, the effective anisotropy γ ef f is manifestly smaller than the intrinsic γ. It increases with temperature, and is in excellent agreement with the experimental data of Fig. 4, strongly suggesting that the modified core cut-off length originating from the tilting of the PV stacks determines the behavior of the melting line at high parallel fields. Very similar behavior has been observed in YBa 2 Cu 3 O 7-δ [16]. Next, the prefactor H ⊥ m0 = B m (0)/µ 0 in Eq. ( 1) is to be interpreted as the hypothetical vortex melting field H m,J (θ = 0) in the absence of the magnetic coupling between PVs. The difference ∆H mag = H ⊥ m (θ = 0) -H m,J (0) ≈ 0.15H ⊥ m between the real (experimental) melting field and this prefactor represents the (remarkably modest) enhancement of the melting field due to magnetic coupling. Summarizing, we have established the existence of phase transition of the vortex lattice in single crystalline Bi 2 Sr 2 CaCu 2 O 8+δ in oblique fields. The transition has a strong first order character and intersects the usual first order vortex lattice melting line at a tricritical point [H ⊥ m (T ), H ct (T )]. For fields parallel to the superconducting layers H < H cr the melting line shows the signature of a composite lattice. For H > H cr , the melting line is fully consistent with that of a uniformly tilted lattice of PV stacks. We thus propose that the new first order transition takes place between the combined and the tilted vortex lattice. For low in-plane fields, the combined vortex lattice is stabilized by the magnetic interaction between PV's in the same stack (vortex line). The enhancement of the melting line in the combined lattice regime is due to the contribution of this magnetic interaction to the vortex line tilt stiffness.

FIG. 1 :

 1 FIG. 1: (a) dc local magnetization loops recorded on an asgrown Bi2Sr2CaCu2O 8+δ single crystal (Tc = 88 K) at 75 K, as function of the magnetic field component H ⊥ perpendicular to the superconducting layers, with the in-plane field H held constant. The inset shows a magnified view of the discontinuity at the vortex melting transition in H = 0. (b) the in-phase (screening) component of the ac response of the same crystal, recorded under the same conditions, with an ac magnetic field of amplitude hac = 0.8 Oe and frequency f = 11 Hz applied along the c-axis. The melting transition shows up as a paramagnetic peak (see inset), the transition from combined to tilted vortex lattice is indicated by arrows. The ac response is plotted as the transmittivity T ′ ≡ [B ′ (f, T ) -B(f, T ≪ Tc)]/[B(f, T > Tc) -B(f, T ≪ Tc)].

FIG. 2 :

 2 FIG.2: Transmittivity T ′ of the same crystal as in Fig.1, recorded at T = 70 K and constant H ⊥ = 58 Oe, as a function of H for various frequencies f and amplitudes hac of the ac ripple field. A marked hysteresis of the ac screening is observed. This hysteresis disappears when hac is increased.

FIG. 3 :

 3 FIG.3: Two field-component vortex-lattice phase diagram, with the first-order melting transition H ⊥ m (•) determined from the paramagnetic peak in T ′ , and the first-order transition to the tilted PV lattice at Bt (•), determined from the "glitch" in T ′ (Fig.1 b). The dashed line is a fit to the composite-to-tilted lattice transition, Eq. (2) with C = 0.030; the continuous line is a fit of the high-field portion of the vortex lattice melting line to Eq. (1).


Wided Miled
email: miled@telecom-paristech.fr

Béatrice Pesquet-Popescu

Jean-Christophe Pesquet
email: pesquet@telecom-paristech.fr

A Convex Programming Approach for Color Stereo Matching

This paper addresses the problem of dense disparity estimation from a pair of color stereo images. Based on a convex set theoretic formulation, the stereo matching problem is cast as a convex programming problem in which a color-based objective function is minimized under specific convex constraints. These constraints arise from prior knowledge and rely on various properties of the disparity field to be estimated. The resulting multi-constrained optimization problem is solved via an efficient parallel block-iterative algorithm. Four different color spaces have been tested in order to evaluate their suitability for stereo matching. Experiments on standard stereo images show that the matching results have been efficiently improved when using color information instead of grey values.

I. INTRODUCTION

Stereo matching is a crucially important problem in computer vision with a wide range of applications, including multiview video coding, multi-viewpoint generation, safe navigation and 3D television. The goal of stereo matching is to compute the disparity map between a pair of stereo images taken from distinct viewpoints. The disparity is defined as the difference in location between corresponding pixels, i.e. pixels resulting from the projection of the same 3D point onto the two image planes.

A broad range of approaches have been developed for solving the stereo matching problem. They are basically classified into three categories : feature-based, area-based and energybased approaches. A good survey for the different stereo matching strategies is addressed in [START_REF] Scharstein | A taxonomy and evaluation of dense twoframe stereo correspondence algorithms[END_REF]. Feature-based methods [START_REF] Veksler | Semi-dense stereo correspondence with dense features[END_REF] provide sparse disparity maps by matching extracted salient features from both images, such as edges, corners or segments. They have a low computational complexity and are well suited for real time processing tasks. In addition, they allow the recovery of large displacements and establish accurate disparity estimates. However, methods computing only sparse matches cannot be considered in many applications of stereo, such as view synthesis and 3D reconstruction. On the other hand, areabased methods produce dense disparity maps by correlation over local windows. They perform well in highly textured regions but often produce noisy disparities in textureless areas and fail around depth discontinuities. Many attempts have been made to solve these problems by varying the size and shape of the matching window according to the local variation of disparity characteristics [3], [START_REF] Fusiello | Symmetric stereo with multiple windowing[END_REF], [START_REF] Davis | Window-based, discontinuity preserving stereo[END_REF]. Energy-based methods provide dense disparity fields by optimizing a global energy function, which is typically the sum of a data term and a smoothness term [START_REF] Bleyer | A layered stereo algorithm using image segmentation and global visibility constraints[END_REF], [START_REF] Slesareva | Optic flow goes stereo: A variational method for estimating discontinuity-preserving dense disparity maps[END_REF], [START_REF] Veksler | Stereo correspondence by dynamic programming on a tree[END_REF], [START_REF] Kolmogorov | Computing visual correspondence with occlusions using graph cuts[END_REF], [START_REF] Miled | Disparity map estimation using a total variation bound[END_REF]. They are called global approaches and are generally more accurate than area-based approaches, especially in the challenging image regions.

In recent years, global optimization approaches have attracted much attention in the stereo vision community due to their excellent experimental results [START_REF] Scharstein | A taxonomy and evaluation of dense twoframe stereo correspondence algorithms[END_REF]. Many global stereo algorithms have, therefore, been developed dealing with ambiguities in stereo such as occlusions, depth discontinuities, lack of texture and photometric variations. These methods exploit various constraints on disparity such as smoothness, visibility, view consistency, etc., while using efficient and powerful optimization algorithms. Although many promising results have been obtained, an even accuracy may be achieved by using the color information, which is typically available in the stereo images.

Color provides much more distinguishable information than intensity values. Therefore, using color images in stereo matching yields more accurate estimates of disparity than grey value images. A number of different approaches have been proposed in the literature for color stereo matching. El Ansari et al. [START_REF] Ansari | A new regions matching for color stereo images[END_REF] proposed a new region based method for matching color images based on the fact that regions contain much richer information than individual pixels. To guarantee a great similarity between corresponding regions, a color-based cost function that takes into account the local properties of region boundaries is used. Koschan et al. [START_REF] Koschan | Color stereo vision using hierarchical block matching and active color illumination[END_REF] used a combination of a hierarchical block matching technique with active color illumination to improve the quality of the matching results, especially in homogenous regions. In [START_REF] Mühlmann | Calculating dense disparity maps from color stereo images, an efficient implementation[END_REF], Mühlmann et al. have developed a real time and efficient method that calculates a correlation based dense disparity map from color stereo images. Alvarez and Sánchez [START_REF] Alvarez | 3-D geometry reconstruction using a color image stereo pair and partial differential equations[END_REF] proposed a generalization of their work presented in [START_REF] Alvarez | Dense disparity map estimation respecting image discontinuities: A PDE and scale-space based approach[END_REF], where they applied a PDEbased method for disparity estimation to color images. They modify the cost function so that it includes all three color components. Other color stereo techniques have been proposed in [START_REF] Koschan | Improving robot vision by color information[END_REF], [START_REF] Belli | Colour contribution for stereo image matching[END_REF], [START_REF] Jordan | Computational stereo vision using color[END_REF] and all have shown that the results have always been improved when using the color information. Most of the methods mentioned above commonly use the RGB color space. However, prior color evaluation studies [START_REF] Koschan | Dense stereo correspondence using polychromatic block matching[END_REF], [START_REF] Chambon | Color stereo matching using correlation measures[END_REF], [START_REF] Bleyer | Evaluation of different methods for using colour information in global stereo matching approaches[END_REF] have shown that the precision of color stereo matching is improved only when a suitable color system is chosen.

In this paper, we propose a global variational method for computing a dense and accurate disparity map from a pair of color images. This method constitutes an extension of the technique we proposed in [START_REF] Miled | Disparity map estimation using a total variation bound[END_REF] to color stereo images. The stereo matching problem is solved through the minimization of a global objective function which is the sum of intensity differences over the three color channels in a predefined color space. Convex constraints introduced in [START_REF] Miled | Disparity map estimation using a total variation bound[END_REF] modelling prior knowledge on the disparity field remain available, except the Nagel-Enkelmann constraint which involves the left stereo image. A generalization of this constraint in the case of color images is proposed, inspired from [START_REF] Alvarez | 3-D geometry reconstruction using a color image stereo pair and partial differential equations[END_REF]. Within a set theoretic framework, the stereo matching problem, formulated as a constrained optimization problem, is then solved via a parallel block-iterative decomposition method. Matching results with different color spaces are finally compared in order to determine the best performing color system.

The rest of the paper is organized as follows. In Section II, we introduce the color stereo model and describe the set theoretic formulation of the disparity estimation problem. Section III is devoted to presenting the objective function we minimize and the constraints we incorporate to the problem. In Section IV, we review the parallel block-iterative algorithm that will be employed to solve this problem. Experimental results on real stereo data sets are illustrated in Section V, followed by a conclusion in Section VI.

II. PROBLEM FORMULATION

A. Color stereo model

The stereo matching problem is to estimate a 2D disparity field by searching for every pixel in the left image the corresponding pixel in the right image. When stereo images are rectified [START_REF] Fusiello | A compact algorithm for rectification of stereo pairs[END_REF], the vertical component of the disparity vector vanishes, so that only a scalar value has to be estimated. Let I l and I r be the left and right color images of a stereo pair, respectively. A color image may be represented as I = I (1) , I (2) , I (3) , where I (k) , k ∈ {1, 2, 3}, represents the k th color channel in the selected color system. Finding a corresponding pixel in the right image I r for each pixel in the left image I l amounts to search the disparity u that minimizes the following cost function, based on the sum of color differences:

J(u) = 3 k=1 (x,y)∈D [I (k) l (x, y) -I (k) r (x -u(x, y), y)] 2 , ( 1
)
where D ⊂ N 2 is the image support. This expression is nonconvex with respect to the displacement field u. Thus, in order to avoid a non-convex minimization, we consider a Taylor expansion of the non-linear terms

I (k) r (x-ū, y), k ∈ {1, 2, 3}
around an initial estimate ū as follows:

I (k) r (x-u, y) ≃ I (k) r (x-ū, y)-(u -ū)∇I (k) r,x (x-ū, y) , (2)
where ∇I (k) r,x (x-ū, y) is the horizontal gradient of the warped right image channel I (k) r . Note that for notation concision, we have not made anymore explicit that u and ū are functions of (x, y) in the above expression.

With the approximation (2), the cost function J under the minimization in (1) becomes quadratic in u. Thus, setting s = (x, y) the spatial position in either image, the objective function to be minimized can be rewritten as:

J(u) = 3 k=1 s∈D [L (k) (s) u(s) -r (k) (s)] 2 (3)
where

L (k) (s) = ∇I (k) r,x (x -ū(s), y) r (k) (s) = I (k) r (x -ū(s), y) + ū(s) L (k) (s) -I (k) l (s).

B. Set theoretic formulation

Minimizing the objective function (3) aims at recovering the best estimate of the disparity image u from the observed fields {L (k) } k and {r (k) } k . This inverse problem is ill-posed due to the fact that the components of {L (k) } k may, simultaneously, vanish. Thus, to convert this problem to a well-posed one, it is useful to incorporate additional constraints modelling prior knowledge and available information on the solution. In the field of computer vision, such constraints were most commonly formulated as additional penalty terms in the objective function. In this work, the problem is addressed from a set theoretic formulation, where each constraint is represented by a convex set in the solution space and the intersection of these sets, the feasibility set, constitutes the family of possible solutions. The aim then is to find an acceptable solution minimizing the given objective function. A formulation of this problem in a Hilbert image space H is therefore:

Find u ∈ S = m i=1 S i such that J(u) = inf J(S) , (4)
where the objective J : H → (-∞, +∞] is a convex function and the constraint sets (S i ) 1≤i≤m are closed convex sets of H. Constraint sets can generally be modelled as level sets:

∀i ∈ {1, . . . , m}, S i = {w ∈ H | f i (w) ≤ δ i } , (5)
where, for all i ∈ {1, . . . , m}, f i : H → R is a continuous convex function and (δ i ) 1≤i≤m are real-valued parameters such that S = m i=1 S i = ∅. The advantage of the convex program (4) is that a wide range of constraints modelling prior information can be explicitly incorporated to the problem as closed convex sets of the form [START_REF] Davis | Window-based, discontinuity preserving stereo[END_REF]. A further advantage of this formulation is to benefit from the availability of powerful optimization algorithms. For the current work, we employ the constrained quadratic minimization method recently developed in [START_REF] Combettes | A block iterative surrogate constraint splitting method for quadratic signal recovery[END_REF] and particularly well adapted to our needs.

III. COLOR STEREO MATCHING

In this section, we introduce the objective function as well as the considered convex constraints to solve the stereo matching problem from color stereo images, within the framework described above.

A. Global objective function

The objective function to be minimized is the quadratic measure (3) derived from our linearized color stereo model (see Section II). To cope with large deviations from the data model, occlusion points which are pixels only visible from one view of the stereo images have been detected and discarded in the expression of the similarity criterion. For a review of various approaches for finding occlusions, we refer to [START_REF] Geoffrey | Detecting binocular half-occlusions: Empirical comparaisons of five approaches[END_REF]. In this work, occlusion points were detected by enforcing the uniqueness and ordering constraints [START_REF] Yuille | A generalized ordering constraint for stereo correspondence[END_REF]. Furthermore, according to the conditions of convergence of the algorithm we use (see Section IV), the objective function J must be strictly convex. However, since the components of {L (k) } k may vanish in (3), J is not secured to be strictly convex. We therefore introduce an additional strictly convex term as follows:

J(u) = 3 k=1 s∈D\O [L (k) (s) u(s)-r (k) (s)] 2 +α s∈D [u(s)-ū(s)] 2 , (6)
where O denotes the occlusion field, ū is an initial estimate and α is a positive constant that weights the second term relatively to the first. The initial estimate ū is obtained, first, from a correlation based method and, then, iteratively refined by choosing the result from a previous estimate as the initial value of the next step. This helps improving the quality of the solution while reducing the sensitivity of the final result to the initial estimate.

B. Convex constraints 1) Disparity range constraint:

The most common constraint on disparity is the knowledge of its range of possible values. Indeed, disparity values are nonnegative and often have known minimal and maximal amplitudes, denoted respectively by u min ≥ 0 and u max . The associated set is

S 1 = {u ∈ H | u min ≤ u ≤ u max } . ( 7
)
2) Total variation regularization constraint: The smoothness constraint, initially introduced by Tikhonov [START_REF] Tikhonov | Solution of ill-posed problems[END_REF], has been one of the most popular regularity assumptions. However, the Tikhonov regularization, by considering a quadratic function, tends to oversmooth discontinuities. In disparity estimation, we are interested in a regularization process that avoids smoothing around object boundaries. This can be achieved with the help of a suitable regularization constraint. In this work, we use a Total Variation based regularization constraint. Initially introduced by Rudin, Osher and Fatemi [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], this regularity measure has emerged as an effective tool to recover smooth images in variational image recovery [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF], which naturally motivates its extension to the field of variational stereo methods [START_REF] Slesareva | Optic flow goes stereo: A variational method for estimating discontinuity-preserving dense disparity maps[END_REF], [START_REF] Miled | Disparity map estimation using a total variation bound[END_REF]. For a differentiable analog image defined on a domain Ω, the total variation is given by: tv

(u) = Ω |∇u(s)| 2 ds , (8)
where ∇u is the gradient of u and | • | 2 denotes the Euclidean norm in R 2 . Practically, tv(u) represents a measure of the lengths of the level lines in the image. In a previous work [START_REF] Miled | Disparity map estimation using a total variation bound[END_REF], we have shown that the total variation tv(u) of the original disparity image u does not exceed some known bound τ . This constraint, which is associated with the set

S 2 = {u ∈ H | tv(u) ≤ τ } , (9)
appears to be particularly relevant in stereo matching, as it smoothes homogenous regions in the disparity image while preserving edges. In addition, the upper bound τ can be estimated with good accuracy from prior experiments and the considered minimization method is shown to be robust with respect to the choice of this bound [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF], [START_REF] Miled | Disparity map estimation using a total variation bound[END_REF].

3) Nagel-Enkelmann constraint: To benefit from the ability of the employed optimization algorithm to incorporate multiple constraints, we use another regularization constraint based on the Nagel-Enkelmann operator [START_REF] Nagel | An investigation of smoothness constraints for the estimation of displacement vector fields from image sequences[END_REF]. The orientedsmoothness ability of this operator has been primarily used in optical flow algorithms [START_REF] Nagel | An investigation of smoothness constraints for the estimation of displacement vector fields from image sequences[END_REF], [START_REF] Alvarez | Reliable estimation of dense optical flow fields with large displacements[END_REF] and has also been used for stereo in the variational framework of [START_REF] Alvarez | Dense disparity map estimation respecting image discontinuities: A PDE and scale-space based approach[END_REF]. A formulation of the associated regularization term as a convex constraint set has been proposed in our previous work [START_REF] Miled | Disparity map estimation using a total variation bound[END_REF]. The main idea of this regularization constraint is that discontinuities in the disparity image are preserved accordingly to the edges of the left image I l . An extension of this constraint, in the case of color images, can be defined as follows:

S 3 = {u ∈ H | Ω ∇u(s) ⊤ D(∇ I l )(s) ∇u(s) ds ≤ δ} , ( 10
)
where δ is a positive constant. The Nagel-Enkelmann operator D(∇ I l ) is given by [START_REF] Alvarez | 3-D geometry reconstruction using a color image stereo pair and partial differential equations[END_REF]:

D(∇ I l ) = 1 |∇ I l | 2 2 + 2γ 2    ∂ I l ∂y -∂ I l ∂x ∂ I l ∂y -∂ I l ∂x ⊤ + γ 2 I    ,
where I denotes the identity matrix, γ is the anisotropic diffusion constant and

∇ I l (s) = ∇I (ks) l (s), where k s = arg max k∈{1,2,3} |∇I (k) l (s)| 2 .
The constraint S 3 has an isotropic behavior within uniform areas (|∇ I l | ≪ γ), but at color edges (|∇ I l | ≫ γ) it introduces an anisotropic smoothing to preserve the discontinuities.

IV. OPTIMIZATION ALGORITHM

The objective of this section is to develop a numerical solution to the problem of color stereo matching which has been formulated as a convex optimization problem. A parallel block iterative algorithm will be employed to efficiently minimize the quadratic objective function [START_REF] Bleyer | A layered stereo algorithm using image segmentation and global visibility constraints[END_REF] over the feasibility set S = ∩ 3 i=1 S i .

A. Subgradient projections

Here, we briefly recall the basic facts on subgradient projections which are necessary for our problem. More details can be found in [START_REF] Combettes | A block iterative surrogate constraint splitting method for quadratic signal recovery[END_REF]. The solution space is a real Hilbert space H, endowed with the standard scalar product . | . and the associated Euclidean norm . . Let S i be the nonempty closed and convex subset of H given by [START_REF] Davis | Window-based, discontinuity preserving stereo[END_REF], where f i is a continuous and convex function. For every u ∈ H, f i possesses at least one subgradient at u, i.e., a vector g i ∈ H such that

∀ z ∈ H, z -u | g i + f i (u) ≤ f i (z) . ( 11
)
The set of all subgradients of f i at u is the subdifferential of f i at u and is denoted by ∂f i (u). If f i is differentiable at u, then ∂f i (u) = {∇f i (w)}. Fix u ∈ H and a subgradient g i ∈ ∂f i (u), the subgradient projection P i u of u onto S i is given by:

P i u =    u - f i (u) -δ i g i 2 g i , if f i (u) > δ i ; u, if f i (u) ≤ δ i . (12)
The proposed algorithm activates the constraints by means of subgradient projections rather than exact projections. The former are much easier to compute than the latter, as they require only the availability of a subgradient (the gradient in the differentiable case). However, when the projection is simple to compute, one can use it as a subgradient projection.

In our case, exact projections onto S 1 is straightforwardly obtained, whereas a subgradient projection onto S 3 can be easily calculated. For the constraint S 2 , the expression of a subgradient projection is given in [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF].

B. Proposed algorithm

We now proceed to the description of the proposed block iterative algorithm to estimate the disparity u.

Algorithm 1: Fix ε ∈ (0, 1/m[ and set n = 0. Set

Φ(s) = 3 k=1 L (k) (s) 2 , Ψ(s) = 3 k=1 L (k) (s) r (k) (s), and compute u 0 as u 0 (s) = Ψ(s)+α ū(s) Φ(s)+α , if s ∈ D \ O, ū (s), otherwise.
Take a nonempty index set I n ⊆ {1, . . . , m}.

For every i ∈ I n , set a i,n = P i,nu n where P i,n is a subgradient projection of u n onto S i as in [START_REF] Koschan | Color stereo vision using hierarchical block matching and active color illumination[END_REF]. Choose weights {ξ i,n } i∈In in (ε, 1] such that

i∈In ξ i,n = 1. Set z n = i∈In ξ i,n a i,n and κ n = i∈In ξ i,n a i,n 2 . If κ n = 0, exit iteration. Otherwise, set • b n = u 0 -u n , • c n such that c n (s) = (Φ(s) + α) b n (s), if s ∈ D \ O, α b n (s), otherwise.
• d n such that

d n (s) = (Φ(s) + α) -1 z n (s), if s ∈ D \ O, α -1 z n (s), otherwise. • λ n = κ n / d n , z n . Set dn = λ n d n , π n = -c n , dn , µ n = b n , c n , ν n = λ n dn , z n and ρ n = µ n ν n -π 2 n . Set u n+1 =      u n + dn , if ρ n = 0, π n ≥ 0; u 0 + (1 + πn νn ) dn , if ρ n > 0, π n ν n ≥ ρ n ; u n + νn ρn (π n b n + µ n dn ), if ρ n > 0, π n ν n < ρ n .
Increment n and go to step .

Theorem 2: Suppose that there exists a positive integer M such that

∀ n ∈ N, n+M-1 k=n I k = {1, . . . , m}, (13)
then every sequence (u n ) n generated by Algorithm 1 converges to the unique solution of (4).

Remarks 3:

Algorithm 1 allows to easily incorporate additional convex constraints if these are available. Its ability to use approximate (subgradient) projections onto the constraint sets makes it possible to handle a wide range of complex convex constraints. In addition, it has been shown in [START_REF] Combettes | A block iterative surrogate constraint splitting method for quadratic signal recovery[END_REF] that, due to its block iterative structure, this algorithm offers a lot of flexibility in terms of parallel implementation. In particular, several processors can be used in parallel to compute the subgradient projections on the different constraint sets (S i ) 1≤i≤m , leading to improved results while reducing the computational time.

V. EXPERIMENTAL RESULTS

In this section, we present results of the proposed method using standard data sets taken from the Middlebury Database. Figure 1 shows the four stereo pairs considered in this work along with their ground truth images. To parameterize our method, the parameter α in equation ( 6) and the anisotropic diffusion constant γ were set to 10 and 1, respectively. Bounds on the constraint sets (S i ) 1≤i≤3 were computed directly from ground truth fields (see Table I) and three cycles of iterations were performed to refine the initial disparity fields as described in Section III-A. The algorithm has been tested with the four color models RGB, LUV, LAB, I 1 I 2 I 3 and the grey level image representation. As ground truth fields are available for the considered stereo pairs, we evaluate the different results quantitatively by computing two error measures: the mean absolute error (MAE) between computed and ground truth fields and the percentage of bad matching pixels (Err) with absolute error larger than one pixel. Following [START_REF] Scharstein | A taxonomy and evaluation of dense twoframe stereo correspondence algorithms[END_REF], we only consider non-occluded pixels when computing the disparity errors. The overall results provided by our method for the different color spaces are shown in Table II(a), where we also indicate the rank of the color spaces according to their MAE and Err errors. As we can see, the precision of the matching has generally been improved when using the color information, except for the Venus stereo pair where color seems to slightly worsen the results. For the other stereo pairs, the mean absolute error was significantly reduced when using the LUV, RGB and I 1 I 2 I 3 color spaces. However, no significant changes in the results have been noticed when using the LAB color space instead of the grey value information. Table II also includes, for comparison purposes, the results from the global optimization algorithm (GC) of Kolmogorov and Zabih [START_REF] Kolmogorov | Computing visual correspondence with occlusions using graph cuts[END_REF] based on graph-cuts. The same observation about the obvious utility of color information in solving the stereo matching problem could be made when comparing the results of the grey value based matching of this method and the color based matching. We can also notice from Table II that the proposed method, compared to the (CG) algorithm, leads always to the best results.

In Figures 2 and3, we show the disparity maps computed by the proposed and the (GC) methods for the three stereo pairs : Teddy, Dolls and Baby. These image pairs are more challenging than the Venus stereo pair, since they have complex scene structures, wide disparity ranges and large occluded regions. When only using grey values, the disparity maps are represented in Figures 23). The results of using the RGB and LUV color spaces are shown in Figures 23(c) and 2-3(d), respectively. As expected, many matching errors are reduced by using the color information. Especially, we notice that the most precise results have been generally obtained by using the LUV space, which seems to be a suitable color space for stereo matching.

VI. CONCLUSION

Stereo matching is a crucially important problem in computer vision. In this paper, we presented a convex programming approach for matching color stereo images. A robust and efficient optimization algorithm was developed to provide a reliable dense disparity map. Within a convex set theoretic framework, this algorithm minimizes a quadratic colorbased objective function subject to some appropriate convex constraints. It can incorporate a wide range of additional constraints and is highly parallelisable due to its block iterative structure. Four different color spaces have been used and compared for the evaluation of matching results from color stereo images. The results indicate that stereo matching can be significantly improved by using a suitable color space. We found that the LUV color space offers the best performances for the considered test image pairs.

Fig. 1 .

 1 Fig. 1. Left images (Top) of the considered stereo pairs and corresponding ground truth images (Down). From left to right: Venus, Teddy, Dolls, Baby.

Fig. 2 .

 2 Fig. 2. Disparity maps from our method applied on Teddy, Dolls and Baby stereo pairs.