D Robbes

C Jorel

E Olivier

L Méchin

S Lebargy

R Bouregba

G Poullain

C Cibert

Electric Field Sensing: what is brought by duality from Flux Gates ?

Ways to sense electric field using PZT capacitors are investigated. Duality rules a priori apply to material physics, linking polarization to magnetization properties, i.e. D(E) -B(H) curves, but no real electrical flux gate or charge gate does exist which would be the dual of common flux gate magnetometer. The Landau -Khalatnikov equations are applied to ferroelectric materials subjected to both a DC and AC pumping signal. They are solved to get model of the PZT capacitors dynamics. They are experimentally investigated using the Sawyer Tower method. A second harmonic detection is a way to extract the DC excitation.

INTRODUCTION

Using theMany measuring principles are currently used to manufacture various magnetic sensors, among which the Magnetic Flux Gates (MFG) are considered as reference sensors for room temperature magnetic field H (or flux density B). Besides, the conjugated physical quantity of H, i.e. the electrical field E, is far from being easily measured, especially if E include a DC component. There are many origins to the quasi lack of E field measuring system. As discussed here, they are one of our motivations. A better understanding of the deep reasons of this unbalance between E and H fields measuring processes, while they are so strongly linked, should then help to get new E field sensors. A methodology based on duality rules applied both to the physical quantities of the used measuring principles as well as to their electronic conditioning units is used. Because search coils, built around magnetic cores, are the basic component of MFG, then search capacitors, built around a ferroelectric (FE) body, should be the basic component of the MFG dual. Such a concept of "Electrical Charges Gate" (ECG) is not similar to the one proposed earlier [START_REF] Kaplan | Duality of the electric covering fieldmill and the flux gate magnetometer[END_REF], that was limited to a functional duality, expressed only at the level of Thévenin-Norton representations. Our contribution starts with a modified concept of the capacitive electrometer. Just based on the dual of the Lenz's law, the measuring principle of search coil magnetometers, we call them search capacitors electrometer. Introducing a FE material Pb(Zr0.6Ti0.4)O3 (PZT) [2] in the previous search capacitor allows us to investigate, both numerically and experimentally, the main features of their signal dynamics when subjected to a voltage having AC and DC components. Understanding such a signals is the necessary step to design ECG using dual analogy with MFG.

The paper is ordered as follow: section II details our motivation and introduces section III, in which we go back to fundamental electric and magnetic induction laws with a balanced formal presentation and finally obtain the search capacitor electrometer. Section IV and V are devoted to PZT capacitors measurements and dynamical modeling. The last section discusses and proposes openings to the present work.

II. MOTIVATIONS

At the very beginning of our work was a questioning in the domain of science learning at the college level: why is it so easy to apply and measure magnetic fields, whereas the same is not true for the electric ones, especially when DC values are involved? That question concerns most of physics teachers in charge of lectures on electricity, electromagnetism and vector's algebra fundamentals, especially when they want to illustrates these fundamentals with experiments. Such a questioning was at the center of [START_REF] Robbes | Autour des équations de Maxwell[END_REF], a paper dedicated to about 15 pedagogical experiments around Maxwell equations. That questioning can also be appreciated after reading [START_REF] Timofeeva | Des lois de l'induction aux théorèmes de Thevenin et de Norton[END_REF]. An equivalent paper devoted to electrical induction laws is still missing : we were unable to find one in literature excepted in the last chapter of [5]. In the framework of an industrial development, a directional, wide bandwidth and wide dynamic range capacitive electrometer [START_REF] Robbes | Physical Quantity Measuring Unit and Device for Measuring a Voltage and an Electric Field[END_REF] was built. Besides science learning many other applications exist : electrostatic discharge (ESD) sources localization, high field, systems control and monitoring, non destructive evaluation (NDE) of insulating and composite materials, non invasive medical instruments (EEG & ECG), electromagnetic interferences (EMI), storms and lightning's activity, etc. Most of these applications would benefit from a low noise instrument working in the DC mode and with an extended bandwidth and dynamic range.

III. BACK TO FUNDAMENTALS OF THE SEARCH CAPACITOR ELECTROMETER

One key missing point when developing electrometers or charge amplifiers is a Norton circuit accounting for the electrical influence of an external E field upon free charges of highly conductive metallic electrodes. The Thévenin equivalent circuit of magnetic capture using a search coil is associated to the mesh equation:

0 = + + + dt d Ri dt di L e ext φ ( 1
)
where φext = Scap Bext is the link to the field measured from the voltage drop e across the coil terminals, once the capture area Scap is known. The circuit current i is set to 0 if a high impedance amplifier is used. The dual equation of Eq (1) is easily written as the node equation of Kirshhoff's circuit laws:

0 = + + + dt dQ Gv dt dv C i ext , (2)
but the physical meaning of the variables in (2), especially Qext is still unclear, even to men of the Art. Chapter IV of [5] proposes a detailed explanation of this meaning, on which we focus here with the help of Fig. 1. It shows a conducting ring together with its physical dual. The former is the basic search coil, with a slot to generate a voltage drop caused by magnetic induction variations. The latter is a pair of conducting discs with a short between them. The electric influence phenomena creates a current -dQext/dt where Qext = Scap Dext is associated to the flux of the classical electrical induction field through a surface Scap. The latter can be defined from a contour line [5]. The dual of the slotted ring looks like a conducting Yo-Yo, the bulky discs of which are two electrical charges reservoirs. The short allows charges exchange between them. Subjected to an external E field drive, the surface charge density adapt itself to produce the induced field that cancels the external one both inside the bulky conducting electrodes and in the thin dielectric that surrounds the short between the electrodes. The amplifier to measure that short-circuit current should be a transimpedance amplifier that exhibits an input admittance much higher than those associated to C and G introduced in Eq (2) so that they can be neglected in the operating frequency range. Once the short is connected to the high admittance input amplifier (using properly guarded cables), that E field capturing device no longer appears as a Yo-Yo, but returns to the basic plan capacitor, justifying the search capacitor's name. Introducing a dielectric material in between the electrode pair, especially those exhibiting non linear properties like PZT's, then benefits from the Yo-Yo picture (except the spatial shift of the short): the externally induced charges exchange between the bulky Yo-Yo electrodes is modulated, thanks to the saturation process that is included in these non linear materials. That picture highlights a mechanism for the proposed Electrical Charge Gate, the dual device of the MFG. The necessary step to get ECG is the accurate modeling of the signal dynamics in PZT based capacitors. A way is proposed in the next section.

IV. PZT CAPACITRORS DYNAMICS

To understand which features of the dynamics are needed, let us remind the measuring principle of MFGs. In Fig. 2a, a typical magnetization curve B(H) of a soft magnetic material is plotted, together with a pumping AC signal. Clearly, the DC value, which is just superimposed to the AC one at the input, is mixed with the AC in the induced signal, the harmonics content of which being then linked to the DC value. This modulation principle is used in MFG to be sensitive to DC fields in magnetic cores coils. That measuring principle is thought to be usable for DC electric field as well, using the polarization curve P(E) (ex. Fig. 2), although the modulation should be modified, due to hysteresis behavior of the P(E) in FE. Lead zirconate titanate or PZTis a perovskite material that has been extensively investigated since [7]. For this work Pb(Zr0.6Ti0.4)O3 thin film capacitors were prepared on Pt/TiO2/SiO2 buffered Si substrates using a multitarget sputtering system. Detailed description of the preparation and characterization of these films are in [8], but they are subjected to fatigue degradations [9]. New fabrication routes are possible using oxide SrRuO3 conductive electrodes to overcome that critical issue for sensors operations [START_REF] Jorel | Epitaxial PZT thin films on YSZ-buffered Si (001) substrates for piezoelectric MEMS or NEMS applications[END_REF]. 250 µm × 250 µm Pt top electrodes and 180 nm thick PZT films were used to study the ferroelectric polarization dynamics and analyze different ways of extracting the DC field from the measurements.The equivalent circuit (Fig. 3) proposed in [START_REF] Sivasubramanian | Physical Kinetics of Ferroelectric Hysteresis[END_REF] and [12] was used as the Landau-Khalatnikov dynamical ferroelectric model: to separate the geometric contribution of the equivalent vacuum capacitor parallel, C0 is used. Computed polarization loops were obtained using MATLAB Simulink, as presented in [START_REF] Didier Robbes | PZT thin films capacitors : a way to sense Electric field ?[END_REF]. In order to advance towards actual sensors circuit design, the modeling was then transferred to SPICE. The phenomenological model developed by Miller [START_REF] Miller | Device modeling of ferroelectric capacitors[END_REF], which uses two hyperbolic tangent branches and an intrinsic multiparameters fitting optimization, was tested but, as it does not use or give any dynamical information, it is not of practical interest for our purpose. The classical Sawyer-Tower test circuit (Fig. 3) [START_REF] Sawyer | Rochelle salt as dielectric[END_REF] was used to measure the D(E) polarization loops. Note that P(E) is plotted as the ε0E contribution is negligible in high permittivity dielectric such as PZT. P represents the total polarization of the ferroelectric film. It includes that part of the polarization related to the polarizability of the material assumed linear, and the ferroelectric part of the polarization related to the switching domains, essentially non linear.The following expressions (3) and ( 4) are used for plotting the hysteresis loops:

PZT out in t (t) V - (t) V E(t) = (3) FE out ST FE ST FE FE A V C A Q A Q P(t) D(t) × = = = = , ( 4
)
where tPZT is the PZT film thickness (180 nm), QFE and QST the charge stored by the FE and ST capacitors and AFE, the FE capacitor area (250 µm × 250 µm). QFE and QST are considered equal since we neglect all the conduction currents. The obtained polarization loops showed the usual behavior with increasing Vin sine wave amplitude Fig. 4 [START_REF] Wouters | ultra-thin PZT films for low-voltage ferroelectric non volatile memories[END_REF]. From the high field part of these curves, we are able to extract the spontaneous polarization and the contribution of the linear permittivity εrLin.

Dynamic modeling, as proposed in [12], has been very easily implemented using the graphical programming language Simulink of the MATLAB environment [START_REF] Didier Robbes | PZT thin films capacitors : a way to sense Electric field ?[END_REF]. An implementation using SPICE is now carried out, producing a FE capacitor dipole. Both dynamic simulations results are shown in Fig. 5, which surrounds the experimental ones. It is worth to note that it is not a mathematical fitting with optimized parameters like in the Miller [START_REF] Miller | Device modeling of ferroelectric capacitors[END_REF] type of modeling. Here, two physical constants are extracted from the experimental loop: the capacitance at E=0 (CE=0 from the slope at E=0) and the resistance RS from the dissipative area of the hysteresis. From there, the two constants of the dimensionless equation [START_REF] Sivasubramanian | Physical Kinetics of Ferroelectric Hysteresis[END_REF] are calculated and implemented in the model. The Curves presented in Fig. 5 are directly computed from the chosen RS and CE=0 values. Nevertheless, these two quantities extracted from the data did not allow a satisfying fit with experimental data. As a comparison, another arbitrary set of these two values is shown Fig. 5. RS and CE=0 were changed from respectively 8.5 kΩ and 6.5 nF to 5.4 kΩ and 3.1 nF. Apart from a slightly better fit of the effective coercive field and loop surface, an a posteriori tailoring of the two parameters does not enhance the match between simulated and measured loops. In fact, that agreement is satisfactory, and, most important, simulated and measured temporal responses are very close to each other.

C ST C 0 C FE R S V in V in V out FE capacitor

V. ELECTRICAL DC FIELD EXTRACTION AND CONCLUSION

As presented earlier on time plots of the voltage Vout (Fig. 3) and on its spectrum [START_REF] Didier Robbes | PZT thin films capacitors : a way to sense Electric field ?[END_REF], the second harmonic component of Vout is similar to that used in MFG. Here we present lock-in measurements of Vout and simulation performed with Vin = 10 VPP at 500 Hz. They are presented Fig 6. in which the same two constants extracted from the experimental data are used for the simulation. The conversion coefficient between the DC input field and the AC voltage output is about 0.7 mV/kVcm -1 for experimental plots and about 1.2 mV/kVcm -1 for the simulated ones using the Landau-Khalatnikov modeling. At this level, we consider this disagreement not so large, associated to the modeled characteristics that are less rounded than the experimental ones. Most of the behavior of the needed parametric signal transfer are clearly obtained, they will help the design of the Electric Charge Gate to be investigated in a next step, both with regards to the input coupling and to other outputs conditioning.

Figure 1 .

 1 Figure 1. Idealized search coil magnetometer and sear capacitor electrometer. The sensors physical outputs is conveyed by a line to the high impedance or high admittance inputs amplifiers.

Figure 2 .

 2 Figure 2. Idealized search coil magnetometer and sear capacitor electrometer. The sensors physical outputs is conveyed by a line to the high impedance or high admittance inputs amplifiers.

Figure 3 .Figure 4 .

 34 Figure 3. Sawyer-Tower circuit including the ferroelectric capacitor (FE) equivalent circuit. C0 represents the vaccum capacitor linked to ε0, Rs the dissipation component and CFE, the ferroelectric contribution. CST=22nF was used as the ST (Sawyer-Tower) capacitor to measure the charge on the FE capacitor CFE

Figure 5 .

 5 Figure 5. Measured and simulated hysteresis loop at 500 Hz and with 10 VPP sin wave excitation. Green simulated curves are obtained with a constants arbitrary shift.

Figure 6 .

 6 Figure 6. 2 nd harmonic amplitude vs. external applied static field E. ϕ is the phase difference with the AC applied Vin.


Julien Fontchastagner

Thierry Lubin

Frédéric Messine

Smaïl Mezani

Efficient Design Using Successive Analytical Subproblems Method: Application to Axial Magnetic Couplings

Keywords: Design, Axial magnetic coupling, Global optimization, Successive analytical subproblems

The authors propose a simple but efficient method for design problems. Instead of using directly a heavy numerical code in an optimization process, we successively define some approximated subproblems thanks to an appropriate analytical model. These subproblems can then be solved exactly by a deterministic global optimization solver and the procedure is updated until the approximation matches a 3D FEM computation. The method is applied to an axial magnetic coupling. Comparative results show the effectiveness of this approach, and real design cases are studied thanks to this method.

I. INTRODUCTION

S INCE around thirty years, it becomes usual to formalize design problems as optimization ones, generally defined as MINLP problems (mixed integer non-linear programs) [START_REF] Fitan | The electromagnetic actuator design problem: a general and rational approach[END_REF]. When the behavior of the considered device is defined by a purely explicit analytical model, such design problems can be solved exactly by a global deterministic algorithm (based on interval arithmetic) as in [START_REF] Fitan | The electromagnetic actuator design problem: a general and rational approach[END_REF] or [START_REF] Fontchastagner | Some Co-axial Magnetic Couplings Designed Using an Analytical Model and an Exact Global Opti mization Code[END_REF].

However some criteria or constraints are often defined thanks to a numerical model and the resulting problem can become hard to solve. The authors present in this paper a simple but efficient technique to solve this kind of design problems when an analytical approximation exists.

II. AXIAL MAGNETIC COUPLINGS

A basic structure of an axial magnetic coupling is considered as shown in Fig. 1. It consists in two rings of permanent magnets, with alternatively north and south poles glued on two flat iron cylinders. The whole geometry is defined with six parameters (five real numbers and one integer) listed in Table . I. One real physical parameter is also added: B r the remanence of the magnets. In this work, we focus on the steady-state behavior of the coupling. The two parts, the driving and the driven ones rotate at a synchronous constant speed. Thus, it is assumed that there is no Foucault's current in the magnets or in the yokes.

A. Analytical Modeling

In [START_REF] Lubin | Simple Analytical Expressions for the Force and Torque of Axial Magnetic Couplings[END_REF], the authors develop, test and use a complete analytical explicit model for such a structure. So it is possible to compute the electromagnetic torque very quickly. This model is based on an analytical resolution of the magnetic vector potential by separating the variables. If only the first harmonic is considered, and after simplifications, the electromagnetic torque is provided by:

Γ ana = 16 B 2 r 3πµ 0 (R 3 2 -R 3 1 ) sin 2 α π 2 sinh 2 (a) sin(pδ) sinh(2(1 + ν)a) (1)
with

a = p 2h R 1 + R 2 and ν = e 2h .
µ 0 is the permeability of vacuum and δ is the relative angular position between the magnets of the lower side with those of the upper side. Then the maximum torque value is obtained when δ = π 2p . Such a model can be really useful in order to determine the influence of given quantities with parametric variations as in [START_REF] Lubin | Simple Analytical Expressions for the Force and Torque of Axial Magnetic Couplings[END_REF] or [START_REF] Shin | Design and Analysis of Axial Permanent Magnet Couplings Based on 3D FEM[END_REF]. Nevertheless the obtained torque value is generally overestimated as shown in section II-C.

B. Numerical Computation of the Electromagnetic Torque

A 3D FEM model is much suitable in order to reach a better accuracy on the torque computation. The strong formulation of the problem corresponding to the scalar magnetic potential U m is given by:

div(µ grad U m -B r ) = 0, (2)
where µ is the permeability in the considered region (air, magnet or iron), B r = (0, 0, ±B r ) in the magnets and (0, 0, 0) in the rest of the domain.

The thickness of iron yokes (h y ) is computed from parameters values thanks to a simple analytical formula (3) based on Ampère's law and flux conservation considerations.

h y = B r B c α π h (2 h + e) R 1 + R 2 2 p (3)
where B c is the average wished value of the induction in the yoke, the value used in the present work is:

B c = 1.5 T.
The problem is solved using the general commercial finite element solver FlexPDE [5]. After the resolution, the induction vector b = (b x , b y , b z ) in the whole domain is deduced from U m . The electromagnetic torque can be computed by applying the Maxwell stress tensor on S z , a circular area in the air gap between the magnets with a radius equal to the domain one.

A smart method to avoid numerical noise is to average this surface integral along the air gap. Then, a good approximation of the torque Γ 3D can be given directly by a volume integral on the whole air cylinder between the magnets V a :

Γ 3D = 1 e e/2 -e/2   Sz     x y 0   × 1 µ 0   b x b z b y b z b 2 z -b 2 2     •u z dx dy   dz = 1 µ 0 e Va (y b x -x b y ) b z dx dy dz (4)

C. Models vs. Measurements

In order to determine the validity of the two previous models a comparison with experimental values is performed. The used coupling is the same prototype described in [START_REF] Lubin | Simple Analytical Expressions for the Force and Torque of Axial Magnetic Couplings[END_REF] (with p = 6 and µ = 500 µ 0 H•m -1 ), see Fig. 2. to the measurements, and should be the reference estimation in a rigorous design process as described below.

III. SUCCESSIVE ANALYTICAL SUBPROBLEMS (SAS)

A. Problem Statement

The considered design problems to solve can always be written under the form (S g ) defined by:

(S g )      min x∈D f obj (x) s.t. h i (x) = 0, i ∈ [[1, n h ]] g i (x) ≤ 0, i ∈ [[1, n g ]] (5)
where the objective function f obj can be the total mass M a for instance, x is the vector of chosen design variables which can be real ones (as for some dimensions) or integers (as for the pole pair number) or even categorical ones (as kinds of materials). D is the search domain defined by the bound values of these variables. The h i and g i functions define a set of n h equality and n g inequality constraints.

Such problems are called MINLP. In the case where all the functions are purely explicit and analytical the so-formulated global optimization problems can be solved exactly by a particular algorithm, called IBBA for Interval Branch and Bound Algorithm. It is based on a Branch and Bound technique where the bounds are computed using interval analysis, see [START_REF] Fitan | The electromagnetic actuator design problem: a general and rational approach[END_REF], [START_REF] Fontchastagner | Some Co-axial Magnetic Couplings Designed Using an Analytical Model and an Exact Global Opti mization Code[END_REF] and [START_REF] Messine | Essays and Surveys in Global Optimization, Chap.10 "A Deterministic Global Optimization Algorithm for Design Problems[END_REF] for details.

Nevertheless in the present case, a numerical computation of the torque (in a black-box form) should be preferred to the analytical expression as seen in the previous section. Previous work has shown the possibility to use some blackbox constraints with IBBA. But it works only with numerical constraints (not the objective) and with relative fast models (around 1 second), as described in [START_REF] Fontchastagner | Design of Electrical Rotating Machines by Associating Deterministic Global Optimization Algorithm With Combinatorial Analytical and Numerical Models[END_REF]. In the present case, the CPU time of a torque computation is greatly higher. This is why a new simple way of using global deterministic optimization with black-box constraint or objective is presented below.

B. SAS Method

The idea is to use both advantages of the two modeling tools: swiftness of the analytical and accuracy of the numerical one. A corrected analytical model is defined by multiplying the initial model by a corrective factor β k . This modified model is used as a surrogate one and it makes it possible to compute the torque rapidly. The associated optimization subproblem (S k ) is solved using IBBA. A comparison with 3D FEM for the found minimum is performed, the corrective factor is updated and the procedure iterated until the difference between the two models is less than a fixed tolerance ε tol .

The method is called SAS (for Successive Analytical Subproblems). In the case where, for instance, the torque defines the first inequality constraint (Γ 3D (x) ≥ Γ fixed ), the corresponding algorithm is given by Algorithm 1.

Algorithm 1 SAS: Successive Analytical Subproblems

ε = ∞; β 1 = 1; k = 0 {Initialisation} while (ε > ε tol ) do k = k + 1
Solve exactly subproblem (S k ) with IBBA:

(S k )        min x∈D f obj (x) s.t. β k • Γ ana (x) ≥ Γ fixed g i (x) ≤ 0, i ∈ [[2, n g ]] h i (x) = 0, i ∈ [[1, n h ]] ∼ x k = argmin(S k ); Γ num = Γ 3D ( ∼ x k ) ε = β k • Γ ana ( ∼ x k ) -Γ num Γ num ; β k+1 = Γ num Γ ana ( ∼ x k ) end while ∼ x = ∼ x k {Found solution}
The interest of such a technique lies in its simplicity. Of course, this technique does not permit to ensure with certainty the global exactness of the found solution. Each ∼ x k is sure to be exact; for the final ∼ x, it depends on the goodness of the analytical model. However, it permits to solve the general problem (5) effectively, as shown thereafter.

IV. NUMERICAL EXPERIMENTS

A. Realistic design cases Some real problems of magnetic coupling design are solved in this section. The chosen variables are:

x = (R 1 , R 2 , h, α, p, σ m ) (6)
where σ m ∈ {1, ..., 5} is a categorical variable representing the kind of magnet material, detailed in Table II. The factor p f , given in this table, is a price factor expressed in unit per kilogram. It means that, for example, the 38M grade is 1.4 times more expensive than the 30M one. The air-gap length e is not considered as a design variable, for all the numerical tests it is set to e = 5 mm. The considered problems are all under the following form: where the objective function f obj can be the volume of magnets V m , the global external volume V g , the total mass M a , or even a cost indicator c i given by the product of the mass of the magnets by the price factor p f . g 2 and g 3 are two geometric shape inequality constraints. h 1 is an equality constraint reflecting the fact that the angular width of the magnet must be an integer: it is a constructability constraint imposed by the chosen manufacturer [START_REF] Hkcm R Engineering E.K | Magnets Segments Online Design[END_REF]. The R 1 , R 2 and h variables have also been imposed to be rounded to the nearest tenth of a millimeter directly in the optimization process. The density of the chosen iron is 7730 kg•m -3 . All the problems have been solved and the obtained results are given in Table III. All the optima satisfy the constraint imposed to the torque (Γ 3D ≥ 100 N•m). Depending on the chosen criterion, several topologies can appear with different pole pair numbers and magnet grades. As expected, the three first solutions use the strongest magnets (σ m = 5). The CPU time are relatively small, and less than four 3D FEM computations are required to reach each solution. An interesting fact: the less expensive device is not obtained with the cheapest magnet grade.

             min x∈D f obj (x) s.t. g 1 : Γ 3D ≥ 100 N•m g 2 : R 2 ≤ 0.9 R 1 , g 3 : R 1 (1 -α) π p ≥ e 2 , h 1 : 180 • α p ∈ IN (7)

B. Comparison with other methodologies

In order to benchmark the efficiency of the developed methodology, a comparison with more classical techniques were performed. Two methodologies are tested:

1) Ga&Num: The 3D FE model is directly associated to a genetic algorithm (GA), with 30 individuals evolving during 50 generations.

2) GA&ANA+SQP&Num: The problem is first solved using only the analytical model with a GA (the torque is computed by 0.8 Γ ana ). The obtained solution is then used as the starting point of a sequential quadratic programming (SQP) dealing only with the real variables.

Both methodologies are based on GA because their use is still widespread [START_REF] Mahmoudi | Design, Analysis, and Prototyping of an Axial-Flux Permanent Magnet Motor Based on Genetic Algorithm and Finite-Element Analysis[END_REF], [START_REF] Ishikawa | Optimization of Rotor Topology in PM Synchronous Motors by Genetic Algorithm Considering Cluster of Materials and Cleaning Procedure[END_REF]. The MATLAB R Global Optimization Toolbox implementations have been chosen. The h1 constraint has been removed to ensure the convergence of these algorithms. The obtained results are given in Table IV. The second methodology provides really good results, but it clearly appears that the SAS method is much faster for the same accuracy.

C. Black-box objective case

The SAS method can also be applied to solve problems with a black-box objective function, e.g. in the present case maximize the torque. The algorithm is nearly the same as Algorithm 1, except that f obj is replaced by -β k • Γ ana and, in this example, we chose a price limit for the first inequality constraint: c i (x) ≤ 1. Results are given in Table V and are compared to the exact solution of the same problem using only Γ ana as the objective function. The same optimum is found which means two things. Firstly, the analytical model is quite good and really gives the right tendencies. Secondly, the SAS method can find the global optimum which, here, is the solution given by IBBA.

D. Multi-objective case

It is also possible to apply the SAS technique sequentially in order to solve multiobjective problems. For instance, Fig. 4 deals with the minimization of the cost c i versus the maximization of the torque Γ 3D . Topological groups appear clearly on the obtained Pareto front (57 3D FEM computations, about 10 h). V. CONCLUSION In this paper an original and simple method, based on successive analytical subproblems, for solving design problems is presented and applied to an axial magnetic coupling. Many concrete problems are solved using different techniques and numerical results show the efficiency and the interest of this new method SAS.

Fig. 1 .

 1 Fig. 1. Geometry of the studied structure (only 2 poles)

Fig. 2 .Fig. 3 .

 23 Fig. 2. Axial magnetic coupling prototype placed on the test benchThe torque values measured on the bench and those computed by the two models are plotted on Fig.3. The abscissa corresponds to the angular displacement in mechanical degrees, and is then equal to 180 π δ. The maximum torque value is reached when this angle is equal to 15 degrees which corresponds to δ = π 2 p as seen in II-A. This maximum value is the chosen dimensioning torque value.It appears clearly that the analytical model gives the right tendencies but overestimates significantly the interesting quantity, e.g. the torque. The 3D FEM model is really much closer

Fig. 4 .

 4 Fig. 4. Pareto front: Cost indicator vs. Electromagnetic torque

TABLE IV SAS

 IV VS. SOME OTHER TECHNIQUES: min Ma with Γ 3D ≥ 100 N•m

		GA&Num *	GA&ANA + SQP&Num *	SAS
	R 1	99.173	93.086	94.1
	R 2	20.253	34.838	28.6
	h	4.429	3.798	3.0
	α	0.417	0.524	0.58
	p	14	15	15
	σm	5	5	5
	Ma (kg)	1.574	1.483	1.483
	Γ 3D (N•m)	100.2	100.0	100.9
	n o of 3D Comp.	1533	123	4
	CPU Time	25 h 30'	1 h 18'	10' 15"

* best run over 10, and without h 1 constraint

TABLE V MAXIMIZATION

 V OF THE TORQUE WITH c i ≤ 1

		SAS	IBBA
		(max Γ 3D	(max Γana)
	R 1 , R 2 , h	100.0, 71.0, 4.9
	α, p, σm	0.794 (11 o ), 13, 2
	Γ 3D (N•m)	108.2	-
	Γana (N•m)		141.3
	c i		1.00
	n o of 3D Comp.	2	0
	CPU Time	4h 20'	2h 13'


Armand Joulin
email: armand.joulin@inria.fr

Bing Sing

Kang
email: sbkang@microsoft.com

Recovering Stereo Pairs from Anaglyphs

An anaglyph is a single image created by selecting complementary colors from a stereo color pair; the user can perceive depth by viewing it through color-filtered glasses. We propose a technique to reconstruct the original color stereo pair given such an anaglyph. We modified SIFT-Flow and use it to initially match the different color channels across the two views. Our technique then iteratively refines the matches, selects the good matches (which defines the "anchor" colors), and propagates the anchor colors. We use a diffusion-based technique for the color propagation, and added a step to suppress unwanted colors. Results on a variety of inputs demonstrate the robustness of our technique. We also extended our method to anaglyph videos by using optic flow between time frames.

Introduction

Arguably, the first 3D motion picture was shown in 1889 by William Friese-Green, who used the separation of colors from the stereo pair to generate color composites. Glasses with appropriate color filters are used to perceive the depth effect. Such images generated by color separation are called anaglyphs. An example is shown in top-left part of Figure 1, where the color separation is red-cyan. Here, anaglyph glasses comprising red-cyan filters are used to perceive depth. There is much legacy anaglyph content available with no original color pairs.

In this paper, we show how we can reliably reconstruct a good approximation of the original color stereo pair given its anaglyph. We assume the left-right color separation is red-cyan, which is the most common; the principles of our technique can be applied to other color separation schemes.

Given an anaglyph, recovering the stereo color pair is equivalent to transferring the intensities in the red channel to the right image and the intensities in the blue and green channels to the left one. Unfortunately, the intensity distributions between the three channels are usually very differ-ent, with mappings of intensities typically being many-tomany. This complicates the use of stereo and standard color transfer methods. Since the input is essentially a stereo pair, we have to handle occlusions as well.

Our technical contributions are as follow:

• A completely automatic technique for reconstructing the stereo pair from an anaglyph (image or video). • Non-trivial extensions to SIFT-Flow [START_REF] Liu | Sift flow: Dense correspondence across different scenes[END_REF] to handle registration across color channels. • Non-trivial extensions to Levin et al.'s [START_REF] Levin | Colorization using optimization[END_REF] colorization technique to handle diffusion of only specific color channels per view, and use of longer range influence. We also added the step of suppressing visually inappropriate colors at disocclusions.

We made the following assumptions:

• There are no new colors in the disoccluded regions.

• The information embedded in the anaglyph is enough to reconstruct full rgb 2 . • The images are roughly rectified (with epipolar deviation of up to ±25 pixels).

As we demonstrate through a variety of results, our technique is robust to inadvertent asymmetric occlusion in one camera (e.g., part of a finger in one view but not in the other) as well as blurry and low-light scenes.

Related work

As far as we know, there are two approaches for generating stereo pairs from anaglyphs. The first is a tool called Deanaglyph3 . However, it requires as input not just the anaglyph, but also one color image. Dietz [START_REF] Dietz | Reprocessing anaglyph images[END_REF] proposed a method to recover the original images from anaglyphs captured using a modified camera with color filters. There is unfortunately insufficient information in the paper on how his technique works; however, his results show that the colors often do not match across the reconstructed stereo pair. Another paper which is closely related is that of Bando et al. [START_REF] Bando | Extracting depth and matte using a color-filtered aperture[END_REF]. They use an RGB color filter on the camera lens aperture to obtain three known shifted views. They propose an algorithm to match the corresponding pixels across the three channels by exploiting the tendency of colors in natural images to form elongated clusters in the RGB space (color lines model of [START_REF] Omer | Color lines: image specific color representation[END_REF]), with the assumption of small disparities (between -5 to 10 pixels only).

We now briefly survey topics relevant to our approach: colorization, extraction of dense correspondences, and color transfer.

Colorization

In our approach, we use a diffusion-based colorization method to propagate the colors of reliable correspondences to other pixels. Colorization methods assume that the color of some input pixels is known and infer the color of the other pixels [START_REF] Levin | Colorization using optimization[END_REF][START_REF] Luan | Natural image colorization[END_REF][START_REF] Yatziv | Image and video colorization[END_REF][START_REF] Li | Scribbleboost: Adding classification to edge-aware interpolation of local image and video adjustments[END_REF][START_REF] Welsh | Transferring color to greyscale images[END_REF]9]. Most colorization methods also assume, explicitly or implicitly, that the greyscale values are known and used as part of the optimization. For anaglyphs, the luminance values are not available.

Our method is based on diffusion of reliably transferred colors using a graph construct (e.g., diffusing reliable red intensities in the left view). Our approach is related to [START_REF] Levin | Colorization using optimization[END_REF]. They assume that the input colors given by a user are approximate and may be refined depending on the image content. In our work, this assumption is not required. Unlike [START_REF] Levin | Colorization using optimization[END_REF], since we work in rgb space directly, we must preserve edges while diffusing the color (especially over shadowed regions). We use longer-range information (for each pixel, within a 21 × 21 centered window) to diffuse colors.

Another popular approach to colorize an image is to use the geodesic of an image to connect a pixel with unknown color to an input pixel [START_REF] Yatziv | Image and video colorization[END_REF]. This type of approach works well in practice but implicitly assumes the number of input pixels is small, which is not true in our case.

Dense correspondences

Our approach relies on finding dense good correspondences between both images. Our approach adapted the techniques of SIFT-flow [START_REF] Liu | Sift flow: Dense correspondence across different scenes[END_REF] and optical flow [13] to fit the requirements of our problem.

There are many methods designed to compute dense correspondences between images. Most stereo matching methods have been proven to be successful when the cameras are radiometrically similar [START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF]. These are not applicable in our case, since we are matching the red channel in one view with cyan in the other. There are also approaches that handle changes in exposure [START_REF] Kim | Visual correspondence using energy minimization and mutual information[END_REF]7] or illumination [START_REF] Heo | Illumination and camera invariant stereo matching[END_REF]. However, such approaches typically assume a bijective mapping between the two intensity maps, which is violated for anaglyphs.

Bando et al. [START_REF] Bando | Extracting depth and matte using a color-filtered aperture[END_REF] use Markov random field to register across color channels. As mentioned earlier, it is important to note that their assumption of local cluster shape similarity based on the color lines model [START_REF] Omer | Color lines: image specific color representation[END_REF] is facilitated by small disparities.

More recently, HaCohen et al. [START_REF] Hacohen | Non-rigid dense correspondence with applications for image enhancement[END_REF] have proposed a color transfer method based on local correspondence between similar images. They look for correspondences between images with different light conditions and taken by different cameras. Unlike graph matching based approaches, their approach finds correspondences with no spatial constraints. In our work, the two images to be reconstructed are highly spatially related; ignoring this information leads to incorrect colorization.

Color transfer

A related subject of interest is color transfer, which is the process of modifying the color of a target image to match the color characteristics of a source image. Most approaches [START_REF] Reinhard | Color transfer between images[END_REF][START_REF] Xiao | Color transfer in correlated color space[END_REF][START_REF] Dong | Fast local color transfer via dominant colors mapping[END_REF]17] assume that at least the target image is in color (or at least greyscale). As a result, they are not directly applicable in our case. These approaches tend to make use of global statistics of the image and thus may be less effective in recovering the true color locally.

Overview of our approach

Our approach to recover the stereo pair from an anaglyph is to iteratively find dense good correspondences between the two images and recolorizing them based on these correspondences. As a post-processing stage, we detect "incompatible" colors in unmatched regions (defined as colors that exist in one image but not the other) and assign them to the nearest known colors. This process prevents new colors from appearing in the occluded regions.

Figure 1 summarizes our process. We use a modified version of SIFT-Flow to obtain an initial set of correspondences between the two images. We then recolorize the images Ground truth SIFT-Flow ASIFT-Flow (ours)

Anaglyph [10] [5] Ours

Figure 3. Comparison between different algorithms in terms of reconstruction. We implemented the algorithm described in [START_REF] Kim | Visual correspondence using energy minimization and mutual information[END_REF] and used the online code of [START_REF] Hacohen | Non-rigid dense correspondence with applications for image enhancement[END_REF].

based on these correspondences. Each channel is then iteratively matched independently using an optical flow method to produce new good correspondences; the colors are updated using the new correspondences.

In the next section, we describe the dense matching algorithms to align the color channels of the two images. We then explain how we retrieve good correspondences from this alignment.

Channel alignment

Given a red-cyan anaglyph, we first need to align the red channel with the blue and green channels. To accomplish this, we use two different graph matching algorithms, namely, modified versions of SIFT-Flow and an optical flow algorithm. We chose SIFT-Flow [START_REF] Liu | Sift flow: Dense correspondence across different scenes[END_REF] because it is capable of locally matching images with different intensity distributions.

Initial rough channel alignment

SIFT-Flow is to compute correspondences between the initially color-separated images. The intuition is that despite brightness variability across the color channels, there is usually a strong correlation in local texture. SIFT is used because it is some form of a texture descriptor and it is inherently invariant to some amount of brightness changes.

We denote by L and R the left and right images respectively. Given pixel p = (x, y) of the left image (resp. right image), we denote by S L (p) (resp. S R (p)) its associated SIFT feature. SIFT-Flow maps pixels of a source image to those in a target image while maintaining the local structure of the input image. More precisely, given the left image as the input, for each pixel p, SIFT-Flow tries to find the best disparity

d(p) = (d x (p), d y (p)) as to match S L (p) to S R (p + d(p)) by minimizing E SF = p∈L S L (p) -S R (p + d(p)) 1 + 1 σ 2 p d(p) 2 + p∼q min(α|d x (p) -d x (q)|, t x ) + min(α|d y (p) -d y (q)|, t y ), (1)
where σ, α are constants and t x and t y are thresholds on the maximum disparity beyond which there is a penalty. We set t x = 5, t y = 2, σ = 14 and α = 2. SIFT-Flow is originally used to match scenes based on context. In our case, the scenes to be matched are the same, just shown as shifted versions and in different color bands. Unfortunately, the original SIFT-Flow tends to not do well in less textured regions, which often results in artifacts if used as is for matching followed by color transfer. These two observations led us to adapt SIFT-Flow to our problem (which we will now refer to as "anaglyph SIFT-Flow" or ASIFT-Flow).

ASIFT-Flow. SIFT flow is designed to match very different images. Each given pixel p is associated with a search window W (p). Typically, in the original SIFT-Flow, this window is square and cover most of the image. In our case, the images are very similar and roughly rectified, and as a result shifts are along the x-axis. We thus restrict the window W (p) to a w 12 × w rectangle along the x-axis, where w is the width of the image (typically w = 640, leading to a 50 × 640 search window).

SIFT-Flow directly matches SIFT features which are known to have a higher value on highly texture regions which encourage matching of highly texture regions over low texture ones. This introduces strong alignment errors in dark or blurry regions (see Figure 2). Notice the artifacts for SIFT-Flow at the darker lower region caused by incorrect color transfer.

We handle this problem by using a mapping φ related to a Hellinger kernel on the SIFT features. More specifically, given a SIFT feature S(p) = (S d (p)) d≤D S of dimension D S , its Hellinger mapping is a D S dimensional vector F (p) with d-th coordinate equal to

F d (p) = S d (p) D S d ′ =1 S d ′ (p) . ( 2
)
We replace in Eq. 1 the SIFT features by their Hellinger mapping. This mapping assigns more weight to regions with smaller gradients, leading to our result in Figure 2. Finally, in the original SIFT-flow method, the cost of matching a pixel to the outside of the target image is independent of its proximity to the left and right borders of the image. In other terms, SIFT-Flow penalize the matching of a pixel p to the outside by a constant cost K. In our case, since we expect our image pair to be roughly rectified, disparity is mostly horizontal. As such, we penalize vertical shifts as well as based on the amount of deviation from the left or right borders. More precisely, given a pixel p at a L 1 distance d(p) to the closest border, we replace the constant cost K by a cost depending on this distance, namely, Kd(p).

Comparison with other matching methods. Why did we use ASIFT-Flow and not other methods? In particular, Kim et al. [START_REF] Kim | Visual correspondence using energy minimization and mutual information[END_REF] propose a stereo matching algorithm which works for images with different intensity distribution. However, they assume that the mapping of intensity levels is oneto-one, which is not true in general between color channels. HaCohen et al. [START_REF] Hacohen | Non-rigid dense correspondence with applications for image enhancement[END_REF] recover the various textures of the image but is not robust enough to deal with high intensity changes resulting in wrong color estimations. We find that these two approaches do not work as well as ASIFT-Flow; a representative example is shown in Figure 3.

SIFT is computed on 16 × 16 patches; its coarseness is why we use it only to obtain a first rough alignment in the initialization. In the next section, we explain how we refine the channel alignment.

Refining channel alignment

Using the matches by ASIFT-Flow, we have a first rough colorization of the two images. In many cases, the results look reasonable. However, artifacts can occur in areas with significant lighting and/or texture changes. In addition, since ASIFT-Flow is block-based, artifacts are usually generated at object boundaries.

Unfortunately, since the color reconstruction is still only approximate, image rectification may only be approximate as well. As a result, we allow for some slop in the vertical direction when refining the search, and we use optical flow. Note that our problem is not that of stereo, but rather color reconstruction. We compute optical flow independently on the three channels, namely, the red channel for the left image and the blue and green channels for the right image. While some of the re-estimated matches are better, errors still occur. In the next section, we describe a criteria to select the good matches.

Finding good matches

Given matches from left to right M L →R and those from right to left M R →L , different criteria can be used to select the more reliable ones. We use three different criteria: stereo consistency, texture consistency, and detection of textureless regions. These criteria are based on comparing an original known channel with its reprojected (warped) version (M L →R and M R →L ).

Stereo consistency (or match reciprocity). This refers to the property that for a pair of correspondences, in each direction, the corresponding pixel is the best choice. In other words, given a pixel i in the left image, we say that it is a good match wrt stereo consistency if i == M R →L M L →R (i) .

Since the imprecision in SIFT matching at the pixel level makes it difficult to satisfy this criterion exactly, we use a soft version where a pixel is allow to match a pixel within a t x × t y (= 5 × 2) window. The modified criteria are

|x(i) -x(M R →L M L →R (i) )| ≤ t x , (3)
|y(i) -y(M R →L M L →R (i) )| ≤ t y . (4)
Texture consistency. A match is good if its reprojected (warped) value is similar to the observed. In this case, we are comparing SIFT features: original observed S(p) and reprojected Ŝ(p). The match is deemed good the distance between the two features is small enough, i.e., |S(p) -Ŝ(p)| < t s . Typically, we fix t s = 2.

Textureless regions. Textureless regions are a special case, since both our rough and refined alignment steps tend to fail there. We define a pixel to be "textureless" if its associated SIFT feature has a norm below a certain threshold (set to 1). A match is deemed good if both pixels are textureless.

The criteria allow us to select good candidate matches. We make the strong assumption that the reconstructed colors (anchor colors) using these matches is correct. Our goal is now to propagate anchor colors to the other pixels, which is the colorization step.

Colorization

There are many approaches for colorization [8,[START_REF] Yatziv | Image and video colorization[END_REF][START_REF] Levin | Colorization using optimization[END_REF][START_REF] Chia | Semantic colorization with internet images[END_REF]. The most popular are either based on diffusion [START_REF] Levin | Colorization using optimization[END_REF] or geodesic distance [START_REF] Yatziv | Image and video colorization[END_REF]. We did not consider methods based on geodesic distance because they are not suitable for cases where the number of anchor pixels (i.e., good matches in our case) is very high. Based on our experiments, they occupy between 80-90% of the image.

In this section, we describe a colorization method inspired by Levin et al. [START_REF] Levin | Colorization using optimization[END_REF]; non-trivial modifications are necessary to adapt to our unique problem formulation.

Diffusion method

In Levin et al. [START_REF] Levin | Colorization using optimization[END_REF], given a greyscale image, a similarity matrix W between nearby pixels is constructed, from which a diagonal matrix D containing the sum of each row of W and a normalized Laplacian matrix L = I -D -1 W are computed and used for color diffusion. Denoting by I(i) the intensity of a given pixel, we use W with entry W ij being zero if the two pixels i and j are not within 10 pixels of each other; otherwise it is

W ij = exp(-γ I(i) -I(j) 2 2 )δ( I(i) -I(j) 1 ≤ λ I ),
where δ is the indicator (membership) function and λ I is a fixed threshold set to 25. The anchor colors (y k , "known") are fixed, with the rest (y u , "unknown") needed to be computed. Let L uu (resp. L uk ) be the sub-block matrix of L with row and column indexed by the unknown pixels (respectively unknown and known pixels). We solve

min yu∈[0,1] 1 2 (y T u L uu y u + 2y u L uk y k ),
which is equivalent to the following linear problem:

L uu y u = L uk y k .
If N p is the number of unknown pixels, solving the linear problem is O(N 2 p ). For most images, N p is small and this closed form solution is efficiently computed. For large images or video sequences, N p can be very big and computing the closed form solution is prohibitive. In those cases, we reduce the complexity of the algorithm by considering N r non-overlapping horizontal bands (we use N r = 4). This is guaranteed to converge since the cost function is strongly convex and the complexity is O(N r N p ).

Differences with Levin et al. [11]

In [START_REF] Levin | Colorization using optimization[END_REF], the user-specified colors may be modified. This is reasonable, since they are specified via scribbles, which are only approximate. Their method can be regarded as a one-step gradient descent with a good initialization. In addition, their method works in YUV space because the grayscale image is known. In our case, the grayscale image is unknown, and as such, this color space is not useable. Finally, our similarity matrix is different. More specifically, for [START_REF] Levin | Colorization using optimization[END_REF], the similarity function is different and only immediate 4-neighbors are connected. This is justified in the YUV space but does not work as well in rgb space.

We iterate between colorizing both images and aligning the channels. At the end of this process, it is possible for disoccluded regions to have colors not found elsewhere (we term such colors as "incompatible" colors). In the next section, we describe a method to detect such colors and recolorize them.

Detecting incompatible colors

Our approach uses the anchor colors to fill in the missing colors. However, objects in disoccluded regions may be incorrectly colored. As shown in Figure 4, for each scene, a border region of the leftmost image is unobserved in the other view, resulting in out-of-place colorization. To detect Before After Ground truth such occurrences and subsequently color correct, we make the assumption that there are no new colors in the occluded regions in order to produce visually plausible outputs.

To recolorize, we first detect regions with incompatible colors before we recolorize them with longer range connections. Let us first consider the right image; the same reasoning is applied to the left image. Given pixel p in the right image with color c p = (r p , g p , b p ), the goal is to find if its estimated red value rp is compatible with the red values observed in the rest of the image given the same blue and green values.

While it seems natural to compute P (r p | g p , b p ) on the good matches, the function will not be precise because the samples are small given the wide range of blue and green values. In addition, it is not clear how to define a threshold over the probability values. Instead, we consider the Nadaraya-Watson kernel regressor to estimate the most probable red value the pixel p should have given its current estimate. The estimator is

r N W p = q∈Q K( rp-rq h ) q∈Q K( rp-rq h ) r q , (5)
where K is a kernel with bandwidth h = 0.1 and Q is the set of good matches with the same blue and green values. We represent 255 2 possible blue and green combinations with 51 2 bins to reduce space requirements. We use 255 bins for the red channel. K is the kernel associated with Gaussian white noise. Note that the computation of mean-shift superpixels and non-local mean are also based on this regressor. Pixel p is deemed to be an incompatible color if |r N W p -rp | is above a certain threshold λ. In our work, we set λ = 0.3 for intensity in the normalized range of 0 to 1. We design our colorization algorithm with a longer range than that for diffusion (101 × 101 instead of 21 × 21), and it helps to substantially reduce coloring artifacts. Two representative results are shown in the middle column in Figure 4.

Results

We ran our experiments on a laptop with a single core (1.6 GHz) and 2 GB of RAM. The code is that of MATLAB with certain functions implemented in C++. It takes on average between 4-10 minutes to process an image (which is first scaled to 640 × 480). If we only apply the rough alignment, the colorization and the incompatible color detection (referred as "MSF+col"), the process takes between 2-5 minutes. Note also that this code can be easily parallelized since most of the process are done on the left and right views independently.

To evaluate the quality of our results, we ran our experiments using stereo images and videos collected over about two years by one of the authors using FujiFilm FinePix REAL 3D W1 camera and Aiptek 3D-HD camcorder. Here we have ground truth, with the anaglyphs synthetically generated from the stereo pairs. Note that both capture devices were deliberately not calibrated for this work. We also show results of our method on images available online (including those recently taken of Mars); there is no ground truth associated with these images.

Comparison with Deanaglyph. This tool was released by 3dtv; it reconstructs the right image given the anaglyph and the left image. Since Deanaglyph assumes the left image is known, we show results only for the right image (our algorithm still computes both images) in Figure 5. Note the visible artifacts generated by Deanaglyph; it also crops out some amount of the border pixels.

Comparison with Bando et al. [START_REF] Bando | Extracting depth and matte using a color-filtered aperture[END_REF]. We ran our algorithm on inputs captured using the camera described in [START_REF] Bando | Extracting depth and matte using a color-filtered aperture[END_REF] (we used their image and code to generate their result). Because the blue and green channels are shifted as well, we matching the red channel to the blue and green independently. We then apply the same colorization algorithm. As shown in Figure 6, we obtained comparable results (though with less apparent edge artifacts) with the same set of parameters as used to reconstruct the other anaglyphs.

Quantitative results on our dataset. We use our datasets ("CityA" with 94 images and "CityB" with 88 images of indoor and outdoor scenes) to evaluate the effectiveness of our algorithm. Here we have ground truth, from which the input anaglyphs are generated and processed. The quality measure we compute is Peak Signal-to-Noise Ratio (PSNR): Given an image I and its reconstruction Î, the PSNR measure is equal to P SN R = 20 log 10 (max(I)) -10 log(M SE), where M SE = I -Î 2 2 . The computed PSNR values of our algorithm are shown in Table 1. We see that our method performs better than the baseline. Note that recolorizing occluded regions does Input [START_REF] Bando | Extracting depth and matte using a color-filtered aperture[END_REF] MSF+Col Ours (final) not improve the PSNR significantly, mostly because the areas occupied by occluded regions are typically small (but obvious if colorized incorrectly).

Qualitative results and robustness. Figure 7 shows a variety of results for anaglyphs found from different online sources, including one from a comic book. These examples show that our method can work on low-quality images, and seems robust to compression artifacts. Figure 8 shows that our algorithm behaves well despite asymmetric occlusion (finger in one view but not the other), blur, or low light.

Failure cases. The most representative failure case is shown in Figure 9. Here, our algorithm fails on part of the flag because the white and red colors produces an almost constant red distribution on the left view. This causes partially incorrect transfer of blue and green to the left view. Our algorithm also does not work well with superimposed layers; other problematic cases include large disocclusions and thin structures.

Extension to video

Our algorithm can be easily extended to process video. To reduce flickering effects, take into account temporal information (in the form of optical flow across successive frames). Unfortunately, processing all the pixels at once is both compute-and memory-intensive.

Deanaglyph

Ours

Ground truth For efficiency reasons, we perform rough channel alignment each frame independently. We impose additional criteria on good correspondences, that of temporal match reciprocity: two pixels connected by the optical flow between two consecutive frames of the left view are matched to pixels of the right view that are also connected by the optical flow. Let O R t →t+1 be the optical flow from the frame t to t+1 on the right view. The additional criteria for a pixel i of the t-th frame of the left view are O R t →t+1 (M L →R (i)) == M L →R (O L t →t+1 (i)) as well as the right view counterpart. Our colorization algorithm is extended to videos by computing similarity matrices within a frame W t and between consecutive frames, W t,t+1 . The similarity matrix within a frame is similar to the one used it the case of images. Given a pixel i and its optical flow O t →t+1 (i), the entries W ij t,t+1 of similarity matrix between consecutive frames are zero if in frame t+1, the two pixels O t →t+1 (i) and j are not within 10 pixels of each other. Otherwise W ij t,t+1 is equal to: exp(-γ I t (i) -I t+1 (j) 2 2 )δ( I t (i) -I t+1 (j) 1 ≤ λ I ).

We define our similarity matrix W over the whole sequence as the block matrix with W t on the diagonal and W t,t+1 on the second diagonal. We define the Laplacian matrix L = I -D -1 W as before. We use a temporal window approach where each frame is colorized given its consecutive frames. This block coordinate descent scheme converges to the global minimum since the cost function is strongly convex and is linear in the number of pixels to colorize. Figure 10 shows results for an anaglyph video, which illustrate the importance of temporal information on preventing flickering effects caused by inconsistent colors across time. If the frames are processed independently, the resulting color inconsistencies tend to be local and minor, but noticeable.

Concluding remarks

We have described a new technique for reconstructing the stereo color pair from either an anaglyph image or video. Current techniques cannot be used as is because of the unique problem of matching across different color channels, transferring only specific color channels, and making sure inappropriate colors do not appear in either reconstructed view. We carefully adapted both SIFT-Flow [START_REF] Liu | Sift flow: Dense correspondence across different scenes[END_REF] and Levin et al.'s [START_REF] Levin | Colorization using optimization[END_REF] colorization technique to ensure good performance, which is evidenced in our results.

Future work includes refining object edges (perhaps using principles of matting) and adding object-level reasoning to overcome cases of ambiguity (as in the case of the redand-white flag in Figure 9).

Figure 1 .

 1 Figure 1. Our pipeline: Given an anaglyph, we iteratively find correspondences between the left and right images and recolorize the two images based on these correspondences.

Figure 2 .

 2 Figure 2. Comparison between the original SIFT-Flow and our version. Notice the incorrect colors transferred for the original SIFT-Flow in the darker lower region.

Figure 4 .

 4 Figure 4. Examples of recolorization of occluded regions with incompatible colors.

Figure 6 .

 6 Figure 6. Comparison with Bando et al. [1] on images taken with their camera. "MSF+Col" refers to using our ASIFT-Flow (rough matching) and colorization. The right column is the result of applying our full algorithm. Notice that the results of our technique have less apparent edge artifacts.

Figure 5 .

 5 Figure 5. Comparison with Deanaglyph (http://www.3dtv.at/Knowhow/DeAnaglyph_en.aspx).

  (a) SIFT-Flow, (b) ASIFT-Flow (coarse matching), (c) ASIFT-Flow (coarse matching) and colorization, (d) Ours without handling incompatible colors, (e) Ours. The best number on each row is shown in bold. Examples from this dataset can be found in the supplementary file.

Figure 10 .

 10 Figure 10. Results for an anaglyph video. Top row: three frames reconstructed independently. The highlighted areas (dotted ellipses) have inconsistent colors across time. Bottom row: results using temporal information. The colors are now more consistent. More examples are shown in the supplementary material.

Table 1 .

 1 Quantitative results in PSNR.

	Dataset # Images View SIFT-Flow ASIFT-Flow ASIFT-Flow + Col Ours -Occ Ours
	CityA	94	left	20.8	21.7	28.2	28.5	28.6
			right	21.2	22.0	29.1	29.3	29.5
	CityB	88	left	20.5	22.3	30.2	30.7	30.7
			right	19.9	21.6	30.9	31.3	31.4

WILLOW project-team, Laboratoire d'Informatique de l'Ecole Normale Supérieure, ENS/INRIA UMR 8548.

Certain types of anaglyphs, such as the "optimized" anaglyph, throws away the red channel; it fakes the red channel of the left image by combining green and blue channels.

[START_REF] Dietz | Reprocessing anaglyph images[END_REF] http://www.3dtv.at/Knowhow/DeAnaglyph_en.aspx

Acknowledgements. The work done by the first author was partially supported by the European Research Council (SIERRA and VIDEOWORLD projects).


Michael Rubinstein

Armand Joulin

Johannes Kopf

Ce Liu

Unsupervised Joint Object Discovery and Segmentation in Internet Images

Car" (a) Images downloaded from the Internet (b) Our automatic segmentation results (c) State-of-the-art co-segmentation results [8] Image search Object discovery Figure 1. Image datasets collected from Internet search vary considerably in their appearance, and typically include many noise images that do not contain the object of interest (a small subset of the car image dataset is shown in (a); the full dataset is available in the accompanying material). Our algorithm automatically discovers and segments out the common object (b). Note how no objects are discovered for noise images in (b). Most previous co-segmentation methods, in contrast, are designed for more homogeneous datasets in which every image contains the object of interest, and, therefore, their performance degrades in the presence of noise (c).

Introduction

We consider the task of jointly segmenting multiple images containing a common object. The goal is to label each pixel in a set of images according to whether or not it belongs to the underlying common object, with no additional information on the images or the object class 1 . Such capa-bility can be useful for automatic generation of large-scale training sets for object detectors/classifiers, data-driven image synthesis, as well as for improving image-to-text relevance and image search.

The task of simultaneously segmenting multiple images is known as co-segmentation, where joint segmentation essentially serves as a means of compensating for the lack of supervisory data, allowing to infer the visual properties of the foreground object even in the absence of a priori information about the object or the images.

While numerous co-segmentation methods have been proposed, they were shown to work well mostly on small datasets, namely MSRC and iCoseg, containing salient and similar objects. In fact, in most of the images in those datasets the foreground can be quite easily separated from the background based on each image alone (i.e. without cosegmentation, see Section 4.1).

However, Internet image collections, such as the ones returned by image search engines for a given user query, are significantly larger and more diverse (Figure 1(a)). Not only do the objects in images downloaded from the Internet exhibit drastically different style, color, texture, shape, pose, size, location and view-point; but such image collections also contain many noise images-images which do not contain the object of interest at all. These challenges, as we demonstrate, pose great difficulties on existing cosegmentation techniques (Figure 1(c)). In particular, most co-segmentation methods assume every image contains the object of interest, and hence are unable to handle dataset noise.

In this paper, we propose a novel correspondence-based object discovery and co-segmentation algorithm that performs well even in the presence of many noise images. Our algorithm automatically discovers the common object among the majority of images and computes a binary object/background label mask for each image. Images that do not contain the common object are naturally handled by returning an empty labeling (Figure 1(b), Figure 2).

Our algorithm is designed based on the assumption that pixels (features) belonging to the common object should be: (a) salient, i.e. dissimilar to other pixels within their image, and (b) sparse, i.e. similar to pixels (features) in other images with respect to smooth transformations between the images. Given an input image dataset, we build a large-scale graphical model connecting similar images, where dense pixel correspondences are used to capture the object's visual variability. These correspondences between images allow us to separate the common object from the background and visual noise.

We performed extensive evaluation of our proposed approach. Our algorithm produces state-of-the-art results on the established MSRC and iCoseg co-segmentation datasets 2 , and provides considerable improvement over previous methods on several new challenging Internet datasets containing rigid and non-rigid object categories. Our Internet datasets, ground truth labels and results are available for the research community for further investigation at http://people.csail.mit.edu/mrub/ObjectDiscovery.

Related work

Object Discovery. Object discovery has been intensively studied in computer vision. In a supervised setup, objects were treated as topics and images as documents, and generative models such as Latent Dirichlet Allocation (LDA) and Hierarchical Pitman-Yor (HPY) have been used to learn the distribution and segmentation of multiple classes simultaneously [24,22]. Winn and Jojic [26] propose a generative model for the distribution of mask, edge and color for visual objects with respect to a smooth deformation field. Although good object recovery results were reported, the model is limited to particular views of an object.

Recently, PageRank [START_REF] Jing | VisualRank: Applying pagerank to large-scale image search[END_REF] was used to discover regions of interest in a bounding box representation [START_REF] Kim | Unsupervised detection of regions of interest using iterative link analysis[END_REF], and selfsimilarities were used to discover a common pattern in several images [START_REF] Bagon | Detecting and sketching the common[END_REF]. Although in these works no generative models were used to learn the distribution of visual objects, reliable matching and saliency are found to be helpful for object discovery. The notions of matching and saliency were also successfully applied by Fakor et al [START_REF] Kim | Distributed cosegmentation via submodular optimization on anisotropic diffusion[END_REF] proposed an optimization which scales up to even larger datasets. Vicente et al. [25] introduced the notion of "objectness" to the co-segmentation framework, showing that requiring the foreground segment to be an object often improves co-segmentation results significantly. All these techniques, however, maintain the strong assumption that the object is present in all of the images, which is not true for Internet image collections.

Other methods were proposed to handle images which might not contain the common object, either implicitly [START_REF] Joulin | Multi-class cosegmentation[END_REF] or explicitly [START_REF] Kim | On multiple foreground cosegmentation[END_REF]. In particular, Kim and Xing [START_REF] Kim | On multiple foreground cosegmentation[END_REF] show promising results given additional user input, but do not show significant improvement in the unsupervised setting. It is clear that in the context of image search and web browsing, user input cannot be used.

Co-segmentation was also explored in weaklysupervised setups with multiple object categories [START_REF] Rubinstein | Annotation propagation in large image databases via dense image correspondence[END_REF][START_REF] Kuettel | Segmentation propagation in imagenet[END_REF]. While image annotations may facilitate object discovery and segmentation, image tags are often noisy, and bounding boxes or class labels are usually unavailable. In this work we show that it is plausible to automatically discover visual objects from the Internet using image search alone.

Object Discovery and Segmentation

Let I = {I 1 , . . . , I N } be the image dataset consisting of N images. Our goal is to compute the binary masks B = {b 1 , . . . , b N }, where for each image I i and pixel x = (x, y), b i (x) = 1 indicates foreground (the common object), and b i (x) = 0 indicates background (not the object) at location x.

Recall our assumption that for an object of interest the foreground pixels should be salient, i.e. dissimilar to other pixels within their image, and sparse, i.e. similar to nearest neighbors (with possible changes in color, size and position). We first define terms which capture these two properties, and then combine them within an optimization framework to solve for the most likely labels for all pixels in the dataset. An overview of the algorithm is shown in Figure 2.

Image saliency. The saliency of a pixel or a region in an image can be defined in numerous ways and exten- Figure 2. One of these things is not like the others. An illustration of joint object discovery and segmentation by our algorithm on two small datasets of five images each. The images are shown at the top row, with two images common to the two datasets -the face and horse images in columns 1 and 2, respectively. Left: when adding to the two common images three images containing horses (columns 3 -5), our algorithm successfully identifies horse as the common object and face as "noise", resulting in the horses being labeled as foreground and the face being labeled as background (bottom row). Right: when adding to the two common images three images containing faces, face is now recognized as common and horse as noise, and the algorithm labels the faces as foreground and the horse as background. For each dataset, the second row shows the saliency maps, colored from black (less salient) to white (more salient); the third row shows the correspondences between images, illustrated by warping the nearest neighbor image to the source image; and the fourth row shows the matching scores based on the correspondences, colored from black (worse matching) to white (better matching).

sive research in computer and human vision has been devoted to this topic. In our experiments, we used an offthe-shelf saliency measure-Cheng et al.'s Contrast-based Saliency [START_REF] Cheng | Global contrast based salient region detection[END_REF]-that produced sufficiently good saliency estimates for our purposes, but our formulation is not limited to a particular saliency measure and others can be used. Briefly, Cheng et al. [START_REF] Cheng | Global contrast based salient region detection[END_REF] define the saliency of a pixel based on its color contrast to other pixels in the image (how different it is from the other pixels). Since high contrast to surrounding regions is usually a stronger evidence for saliency of a region than high contrast to far away regions, they weigh the contrast by the spatial distances in the image.

Given a saliency map, M i , for each image I i , we first compute the dataset-wide normalized saliency, M i (with values in [0, 1]), and define the term

Φ i saliency (x) = -log Mi(x). ( 1
)
This term will encourage more (resp. less) salient pixels to be labeled foreground (resp. background) later on.

Pixel Correspondence. To exploit the dataset structure and similarity between image regions, we need to establish reliable correspondences between pixels in different images. This enables us to determine a pixel as background even when it may be very salient within its own image. We do this using SIFT flow [START_REF] Liu | SIFT flow: Dense correspondence across scenes and its applications[END_REF], which has been successfully applied in the past for label propagation [START_REF] Liu | Nonparametric scene parsing via label transfer[END_REF][START_REF] Rubinstein | Annotation propagation in large image databases via dense image correspondence[END_REF]. However, instead of establishing the correspondence between all pixels in a pair of images, as done by previous work, we solve and update the correspondences based on our estimation of the foreground regions. This helps in ignoring background clutter and ultimately improves the correspondence between foreground pixels (Figure 3). Formally, let w ij denote the flow field from image I i to image I j . Given the binary masks b i , b j , the SIFT flow objective function becomes

E (wij; bi, bj) = x∈Λ i bi(x) bj(x + wij(x)) Si(x) -Sj(x + wij(x)) 1 + (1 -bj(x + wij(x))C0 + y∈N i x α w(x) -w(y) 2 , (2)
where Si are the dense SIFT descriptors of image Ii, x → x p is the Lp distance for p = 1 and 2, Λi is image Ii's lattice, N i

x is the neighborhood of x, α weighs the smoothness term, and C0 is a large constant. We then denote by W the set of all pixel correspondences in the dataset: W = ∪ N i=1 ∪I j ∈N i wij. The difference between this objective function and the original SIFT flow [START_REF] Liu | SIFT flow: Dense correspondence across scenes and its applications[END_REF] is that it encourages matching foreground pixels in image Ii with foreground pixels in image Ij. We also use an L2norm for the smoothness term instead of the truncated L1-norm in the original formulation [START_REF] Liu | SIFT flow: Dense correspondence across scenes and its applications[END_REF] to make the flow more rigid in order to surface mismatches between the images. Figure 3(a) shows the contribution of this modification for establishing reliable correspondences between similar images.

For small datasets, we can estimate the correspondences between any pair of images, however for large datasets such computation is clearly prohibitive. Therefore, we first find for each image Ii a set of similar images, Ni, based on global image statistics that are more efficient to compute, and estimate pixel correspondences with those images only. For each image Ii, we fix the size of Ni to the same constant, K. We use the Gist descriptor [START_REF] Oliva | Modeling the shape of the scene: A holistic representation of the spatial envelope[END_REF] in our implementation, and similarly modify it to account for the foreground estimates by giving lower weight in the descriptor to pixels labeled as background. Figure 3(b-c) demonstrate that better sorting of the images is achieved when using this weighted Gist descriptor, which in turn improves the set of images with which pixel correspondences are computed.

Based on the computed correspondences, we define the matching term

Φ i match (x) = 1 |Ni| j∈N i Si(x) -Si(x + wij(x) 1 , ( 3
)
where smaller values indicate higher similarity to the corresponding pixels. Similarly to the saliency, we compute a dataset-wide normalized term (with values in [0, 1]), Φ i match .

Foreground Likelihood. We use the above saliency and matching terms to define the likelihood of a pixel label:

Φ i (x) = Φ i saliency (x) + λ match Φ i match (x), bi(x) = 1, β, bi(x) = 0, ( 4
)
where β is a constant parameter for adjusting the likelihood of background pixels. Decreasing β makes every pixel more likely to belong to the background, thus producing a more conservative estimation of the foreground.

Regularization. We would like the masks bi to be spatially consistent within each image, i.e. neighboring pixels are encouraged to have the same label, subject to the image structures. We thus define the intra-image compatibility between adjacent pixels x, y in image Ii as [START_REF] Liu | Nonparametric scene parsing via label transfer[END_REF]19]

Ψ i int (x, y) = [bi(x) = bi(y)] exp -Ii(x) -Ii(y) 2 2 , ( 5
)
where the indicator function [•] is 1 when its argument is true, and 0 otherwise. We would also like the labeling to be consistent between images, and so we add a term accounting for the inter-image compatibility between a pixel x in image Ii and its corresponding pixel y = x + wij(x) in image Ij:

Ψ ij ext (x, y) = [bi(x) = bj(y)] exp -Si(x) -Sj(y) 1 . (6)
Notice that SIFT features are used for the inter-image similarity metric in Equation 6 whereas RGB intensities are used for the intra-image similarity in Equation 5.

Finally, once we have an estimate of bi, we can learn the color histograms of the background and foreground regions of image Ii, denoted h 0 i and h 1 i , respectively. We also denote hi = h 0 i , h 1 i , and H = ∪ N i=1 hi. We add the term Φ i color (x) accounting for the contribution of the pixel x to the foreground or background color model based on the segmentation estimate bi(x):

Φ i color (x, hi) = -log h b i (x) i (x). (7)
We use 3D histograms in color space (with 64 bins in each dimension) to model the color distributions instead of the Gaussian mixture models used in [19]. By combining all the aforementioned terms, we obtain a cost function, E(B; W, H), for the segmentations B given the correspondences W and the color models H:

E(B; W, H) = N i=1 x∈Λ i Φ i (x) + λ color Φ i color (x, hi) + y∈N i x λ int Ψ i int (x, y) + j∈N i λ ext Ψ ij ext (x, x + wij(x)) . (8)
Optimization. Our algorithm alternates between optimizing the correspondences W (Equation 2), and the binary masks B (Equation 8). Instead of optimizing Equation 8 jointly over all the dataset images, we use coordinate descent that already produces good results. More specifically, at each step we optimize for a single image by fixing the segmentation masks for the rest of the images. Note that our cost function is non-convex and is not guaranteed to reach the global minimum. After propagating labels from other images, we optimize each image using a Grabcutlike [18] alternation between optimizing Equation 8 and estimating the color models hi. The algorithm then recomputes neighboring images and pixel correspondences based on the current foreground estimates, and the process is repeated for a few iterations until convergence (we typically used 5 -10 iterations). , measured as the ratio of correctly labeled pixels (both foreground and background), and compared to state-of-the-art co-segmentation methods (we performed a separate comparison with Object Cosegmentation [25]; see the text and Table 1). Each plot shows the average per-class precision on the left, followed by a breakdown of the precision for each class in the dataset.

Results

We conducted extensive experiments to verify our approach, both on standard co-segmentation datasets and image collections downloaded from the Internet. We tuned the algorithm's parameters manually on a small subset of the Internet images, and vary β to control the performance. Unless mentioned otherwise, we used the following parameter settings: λ match = 4, λint = 15, λext = 1, λ color = 2, α = 2, K = 16. Our implementation of the algorithm is comprised of distributed Matlab and C++ code, which we ran on a small cluster with 36 cores.

We present both qualitative and quantitative results, as well as comparisons with state-of-the-art co-segmentation methods on both types of datasets. Quantitative evaluation is performed against manual foreground-background segmentations that are considered as "ground truth". We use two performance metrics: precision, P (the ratio of correctly labeled pixels, both foreground and background), and Jaccard similarity, J (the intersection over union of the result and ground truth segmentations). Both measures are commonly used for evaluation in image segmentation research. We show a sample of the results and comparisons in the paper, and refer the interested reader to many more results that we provide in the supplementary material.

Results on Co-segmentation datasets

We report results for the MSRC dataset [START_REF] Shotton | Textonboost: Joint appearance, shape and context modeling for multi-class object recognition and segmentation[END_REF] (14 object classes; about 30 images per class) and iCoseg dataset [2] (30 classes; varying number of images per class), which have been widely used by previous work to evaluate co-segmentation performance. Both datasets include human-given segmentations that are used for the quantitative evaluation.

We ran our method on these datasets both with and without the inter-image components in our objective function (i.e. when using the parameters above, and when setting λ match = λext = 0, respectively), where the latter effectively reduces the method to segmenting every image independently using its saliency map and spatial regularization (combined in a Grabcut-style iterative optimization). Interestingly, we noticed that using the inter-image terms had negligible effect on the results for these datasets. Moreover, this simple algorithm-an off-the-shelf, low-level saliency measure combined with spatial regularization-which does not use co-segmentation, is sufficient to produce accurate results (and outperforms recent techniques; see below) on the standard cosegmentation datasets!

The reason is twofold: (a) all images in each visual category in those datasets contain the object of interest, and (b) for most of the images the foreground is quite easily separated from the background based on its relative saliency alone. A similar observation was recently made by Vicente et al. [25], who noticed that their single image classifier outperformed recent co-segmentation methods on these datasets, a finding that is reinforced by our experiments. We thus report the results when the inter-image components are disabled. Representative results from a sample of the classes of each dataset are shown in Figure 5.

Comparison with co-segmentation methods. We compare our results with three previously proposed methods [8, [START_REF] Joulin | Multi-class cosegmentation[END_REF][START_REF] Kim | Distributed cosegmentation via submodular optimization on anisotropic diffusion[END_REF]. For all three methods we used the original implementations by the authors that are publicly available, and verified we are able to reproduce the results reported in their papers when running their code. The per-class precision is shown in Figure 4 and the Jaccard similarities are available in the supplemental material. Our overall precision (87.66% MSRC, 89.84% iCoseg) shows significant [25] on MSRCV and iCoseg. P and J denote the average precision and Jaccard similarity, respectively. The per-class performance and visual results are available in the supplementary material.

improvement over [START_REF] Joulin | Multi-class cosegmentation[END_REF] (73.61% MSRC, 70.21% iCoseg) and [START_REF] Kim | Distributed cosegmentation via submodular optimization on anisotropic diffusion[END_REF] (54.65% MSRC, 70.41% iCoseg).

Comparison with Object Cosegmentation [25]. Vicente et al.'s method [25] is currently considered state-of-the-art on these datasets2 . Their code is not publicly available, however they provided us with the segmentation masks for the subsets of MSRC and iCoseg they used in their paper. We performed a separate comparison with their method using only the subset of images they used. Our method outperforms theirs on all classes in MSRC and 9/16 of the classes in iCoseg (see supplementary material), and our average precision and Jaccard similarity are slightly better than theirs (Table 1). We note that despite the incremental improvement over their method on these datasets, our results in this case were produced by segmenting each image separately using generic, low-level image cues, while their method segments the images jointly and requires training.

Results on Internet Datasets

Using the Bing API, we automatically downloaded images for three queries with query expansion through Wikipedia: car (4, 347 images), horse (6, 381 images), and airplane (4, 542 images). With K = 16 nearest neighbors, it took 10 hours on average for the algorithm to process each dataset.

Some discovery results are shown in Figure 6. Overall, our algorithm is able to discover visual objects despite large variation in style, color, texture, pose, scale, position, and viewing-angle. For the objects under a uniform background or with distinct colors, our method is able to output nearly perfect segmentation. Many objects are not very distinctive from the background in terms of color, but they were still successfully discovered due to good correspondences to other images. For car, some car parts are occasionally missing as they may be less salient within their image or not well aligned to other images. Similarly, for horse, the body of horses gets consistently discovered but sometimes legs are missing. More flexible transforms might be needed for establishing correspondences between horses. For airplane, saliency plays a more important role as the uniform skies always match best regardless of the transform. However the algorithm manages to correctly segment out airplanes even when they are less salient, and identifies noise images, such as that of plane cabins and jet engines, as background, since those have an overall worse matching to other images in the dataset.

For qualitative evaluation, we collected partial human labels for each dataset using the LabelMe annotation toolbox [START_REF] Russell | Labelme: a database and web-based tool for image annotation[END_REF] and a combination of volunteers and Mechanical Turk workers, resulting in 1, 306 car, 879 horse, and 561 airplane images labeled. All labels were manually inspected and refined. In Table 2 we show the precision and Jaccard similarity of our method on each dataset, with and without using image correspondences. The performance on airplane is slightly better than horse and car as in many of the images the airplane can be easily segmented out from the uniform sky background. Image correspondences helped the most on the car dataset (+11% precision, +17% Jaccard similarity), probably because in many of the images the cars are not that salient, while they can be matched reliably to similar car images to be segmented correctly.

Comparison with co-segmentation methods. We also compared our results with the same three state-of-the-art cosegmentation methods as in Section 4.1 by running them on our datasets. Since the competing methods do not scale to large datasets, we randomly selected 100 of the images with available ground truth labels from each dataset. We re-ran our method on these smaller datasets for a fair comparison. We also compared to two baselines, one where all the pixels are classified as background ("Baseline 1"), and one where all pixels are classified as foreground ("Baseline 2"). Table 3 summarizes this comparison, showing again that our method produces much better results according to both performance metrics (ours results are not exactly the same as in Table 2 bottom row, since only subsets of the full datasets are used here). The largest gain in precision by our method is on the airplane dataset, which has the highest noise level of these three datasets. Some visual comparisons are shown in Figure 7 and more are available in the supplementary material.

Discussion and Limitations

Some failures of the algorithm are shown in Figure 6 (last row of each dataset) and the supplemental material. False positives include a motorcycle and a headlight in the car dataset, and a tree in the horse dataset. This indicates that although matching image structures often leads to object-level correspondence, exceptions occur especially when context is not taken into account.

The algorithm also fails occasionally to discover objects with unique views or background. This is because Gist is a global image descriptor, and unique view and background make it difficult to retrieve similar objects in the dataset. Finally, our algorithm makes the implicit assumption of nonstructured dataset noise. That is, repeating visual patterns are assumed to be part of some "common" object. For example, had a dataset of 100 car images contained 80 images of cars and 20 images of car wheels, then using K = 16 neighbor images by our algorithm may result in intra-group connections, relating images of cars to other images of cars and images of wheels with others alike. In such case the algorithm may not be able to infer that one category is more common than the other, and both cars and wheels would be segmented as foreground. Fortunately, the fixed setting of K we used seems to perform well in practice, however in the general case K needs to be set according to what the user considers as "common".

Conclusion

We explored automatic visual object discovery and segmentation from the Internet using one query of an object category. Image datasets resulting from such queries are significantly more diverse and noisy than the ones used to develop and evaluate previous co-segmentation work. The common object often differs drastically in appearance, and a significant portion of the images may not contain the object at all. We demonstrated that existing co-segmentation algorithms do not perform well in such cases, and presented a new algorithm that is able to naturally handle the visual variation and noise in Internet images. We model the sparsity and saliency properties of the common object, and construct a large-scale graphical model to jointly infer a binary mask for each image. We demonstrated improvement over existing cosegmentation techniques on standard co-segmentation datasets and several challenging Internet datasets.

  Comparison between standard and weighted Sift flow. (b) Nearest neighbor ordering (left to right) for the source image in (a), computed with the standard Gist descriptor. (c) Nearest neighbor ordering (bottom row; left to right) for the source image in (a), computed with a weighted Gist descriptor using the foreground estimates (top row).

Figure 3 .

 3 Figure 3. Weighted Gist and Sift flow for improved image correspondence. We use the foreground mask estimates to remove background clutter when computing correspondences (a), and to improve the retrieval of neighbor images (compared to (b), the ordering in (c) places right-facing horses first, followed by leftfacing horses, with the (noise) image of a person last).

  ho us e pl an e sh ee p si gn tr ee al lo on s K en do K en do 2 L iv er po ol M on ks S ta tu eo fL ib er ty T ra ck an dF ie ld W in dm il l W om an S oc ce r W om an S oc ce r2 ba se ba ll be ar 2 br ow n_ be ar ch ee ta h el ep ha nt fe rr ar i go os e gy m na st ic 1 gy m na st ic 2 gy m na st ic 3 he li co pt er pa nd a1 pa nd a2 py ra m id sk at e sk at e2 sk at e3 st on eh en ge ta j_ m ah al Joulin et al. [7] Joulin et al. [8] Kim et al. [11] Ours

Figure 4 .

 4 Figure4. Segmentation accuracy on MSRC (left) and iCoseg (right), measured as the ratio of correctly labeled pixels (both foreground and background), and compared to state-of-the-art co-segmentation methods (we performed a separate comparison with Object Cosegmentation[25]; see the text and Table1). Each plot shows the average per-class precision on the left, followed by a breakdown of the precision for each class in the dataset.

Figure 5 .

 5 Figure 5. Sample results on MSRC (top two rows) and iCoseg (bottom two rows). For each image we show a pair of the original (left) and our segmentation result (right). More results and qualitative comparisons with state-of-the-art are available in the supplementary material.

Figure 6 . 7 .

 67 Figure6. Automatic discovery of cars, horses and airplanes downloaded from the Internet, containing 4, 347, 6, 381 and 4, 542 images, respectively. For each image, we show a pair of the original (left) and the segmentation result (right). Notice how images that do not contain the object are labeled as background. The last row of each dataset shows some failure cases where no object was discovered or where the discovery is wrong or incomplete. Quantitative results are available in Table2, and more visual results can be found in the supplementary material.

  . [5], a work done in parallel to ours, for unsupervised discovery of image categories. Co-segmentation. Co-segmentation was first introduced by Rother et al. [19], who used histogram matching to simultaneously segment the same object in two different images. Since then, numerous methods were proposed to im-

	prove and refine the co-segmentation [16, 6, 2, 8], many of
	which work in the context of a pair of images with the exact
	same object [19, 16, 6] or require some form of user inter-
	action [2, 4].
	These techniques were later extended in various ways.
	Joulin et al. [8] used a discriminative clustering frame-
	work that can handle multiple images, and Kim et al.

Table 1 .

 1 Comparison with Object Cosegmentation

	Method	MSRC P	J	P	iCoseg	J
	Vicente et al. [25]	90.2	70.6	85.34	62.04
	Ours	92.16	74.7	89.6	67.63

Table 2 .

 2 Segmentation accuracy on the Internet datasets, with and without utilizing image correspondences. Next to the name of each dataset is its percentage of noisy images (images that do not contain the object). P denotes precision and J denotes Jaccard similarity. Qualitative results for these datasets are shown in Figure6and the supplementary material.

	Method	Car (7.5%)	Horse (7.8%)	Airplane (16%)
		P	J	P	J	P	J
	Without corr. 72.25 46.10 74.88 50.06 80.53	51.18
	With corr.	83.38 63.36 83.69 53.89 86.14	55.62
	Method	Car (11%)	Horse (7%)	Airplane (18%)
		P	J	P	J	P	J
	Baseline 1	68.91	0	81.54	0	87.48	0
	Baseline 2	31.09 34.93 18.46 19.85 12.52	15.26
	Joulin et al. [8]	58.7	37.15 63.84 30.16 49.25	15.36
	Joulin et al. [9]	59.2	35.15 64.22 29.53 47.48	11.72
	Kim et al. [12] 68.85	0.04	75.12	6.43	80.2	7.9
	Ours	85.38 64.42 82.81 51.65 88.04	55.81

Table 3 .

 3 Comparison with previous co-segmentation methods on the Internet datasets.

We note that while we call our method "unsupervised", we do assume that the input image dataset contains a common visual category. We use "unsupervised" to emphasize that, other than this assumption, the algorithm makes no further use of a priori information such as the common object's class or image annotations.

While writing this paper, Kuettel et al.[START_REF] Kuettel | Segmentation propagation in imagenet[END_REF] managed to improve the state-of-the-art precision on the iCoseg dataset (91.4%).

Acknowledgments. We thank Antonio Torralba for his help in collecting human foreground-background segmentations for our Internet datasets. This work was done while Michael Rubinstein and Armand Joulin were interns at Microsoft Research Redmond. Michael Rubinstein is supported by the Microsoft Research PhD Fellowship. Armand Joulin is supported by the European Research Council (SIERRA and VIDEOWORLD projects).