Éric Grégoire

Bertrand Mazure

Cédric Piette

Using Local Search to find MSSes and MUSes

Keywords: SAT, MSSes, MUSes, satisfiability, hybrid algorithm

In this paper, a new complete technique to compute Maximal Satisfiable Subsets (MSSes) and Minimally Unsatisfiable Subformulas (MUSes) of sets of Boolean clauses is introduced. The approach improves the currently most efficient complete technique in several ways. It makes use of the powerful concept of critical clause and of a computationally inexpensive local search oracle to boost an exhaustive algorithm proposed by Liffiton and Sakallah. These features can allow exponential efficiency gains to be obtained. Accordingly, experimental studies show that this new approach outperforms the best current existing exhaustive ones.

Introduction

This last decade, the SAT problem, namely the issue of checking whether a set of Boolean clauses is satisfiable or not, has received much attention from the AI research community. Indeed, SAT appears to be a cornerstone in many domains, like e.g. nonmonotonic reasoning, automated reasoning, model-based diagnosis, planning and knowledge bases verification and validation. However, only knowing that a SAT instance is unsatisfiable is often not satisfactory since we might prefer knowing what goes wrong with the instance when this latter one is expected to be satisfiable.

In this respect, the MUS (Minimally Unsatisfiable Subformula) concept can be crucial since a MUS can be seen as an irreducible cause for infeasibility. Indeed, a Email addresses: gregoire@cril.fr ( Éric Grégoire), mazure@cril.fr (Bertrand Mazure), piette@cril.fr (Cédric Piette).

MUS is an unsatisfiable set of clauses such that any of its subsets is satisfiable. It thus provides one explanation for unsatisfiability that cannot be made shorter in terms of the number of involved clauses. Restoring the satisfiability of an instance cannot be done without fixing all its MUSes.

Unfortunately, a same instance can exhibit several MUSes. Actually, the number of these MUSes can be exponential since a n-clauses SAT instance can exhibit C n/2 n MUSes in the worst case. Moreover, computing MUSes is intractable in the general case. Indeed, checking whether a set of clauses is a MUS or not is DPcomplete [START_REF] Papadimitriou | The complexity of facets resolved[END_REF] and checking whether a formula belongs to the set (clutter) of MUSes of an unsatisfiable SAT instance or not, is in Σ p 2 [START_REF] Eiter | On the complexity of propositional knowledge base revision, updates and counterfactual[END_REF]. Fortunately, the number of MUSes remains often tractable in real-life applications. For example, in modelbased diagnosis [START_REF]Readings in Model-Based Diagnosis[END_REF], it is often assumed that single faults occur most often, which can entail small numbers of MUSes.

A dual concept is the notion of Maximal Satisfiable Subset (MSS) of a SAT instance, and the complement of a MSS in a SAT instance is called a CoMSS. The complete set of MUSes or MSSes is an implicit encoding of the other [START_REF] Liffiton | On finding all minimally unsatisfiable subformulas[END_REF]. Specifically, a CoMSS is a hitting set of the set of MUSes and represent minimal sets of clauses that should be dropped in order to restore consistency. In this paper, we are interested in exhaustive approaches to compute these three correlated concepts in the full Boolean clausal framework.

Recently, several approaches have been proposed to approximate or compute MUSes and MSSes, both in the Boolean framework and for other types of constraints. Some of them concern specific classes of clauses or remain tractable for small instances, only. Among them, let us mention the approach in [START_REF] Bruni | On exact selection of minimally unsatisfiable subformulae[END_REF], where it is shown how a MUS can be extracted in polynomial time through linear programming techniques for clauses exhibiting a so-called integral property. However, only restrictive classes of clauses obey such a property (mainly Horn, renamable Horn, extended Horn, balanced and matched ones). Let us also mention [START_REF] Büning | On subclasses of minimal unsatisfiable formulas[END_REF][7] [START_REF] Fleischner | Polynomial-time recognition of minimal unsatisfiable formulas with fixed clause-variable difference[END_REF], which contain studies of the complexity and the algorithmic aspects of extracting MUSes for specific classes of clauses. In [START_REF] Bruni | Approximating minimal unsatisfiable subformulae by means of adaptive core search[END_REF], an approach is proposed that approximates MUSes by means of an adaptative search guided by clauses hardness. In [START_REF] Zhang | Extracting small unsatisfiable cores from unsatisfiable boolean formula[END_REF] a technique is described, that extracts MUSes by learning nogoods involved in the derivation of the empty clause by resolution. In [START_REF] Lynce | On computing minimum unsatisfiable cores[END_REF], a complete and exhaustive technique to extract smallest MUSes is introduced. In [START_REF] Oh | AMUSE: a minimallyunsatisfiable subformula extractor[END_REF], a DPLL-oriented approach has been presented that is based on a marked clauses concept to allow one to approximate MUSes. In [START_REF] Grégoire | Extracting MUS[END_REF], Grégoire, Mazure and Piette have proposed a heuristic-based incomplete approach to compute MUSes, which outperforms competing ones from a computational point of view. Interestingly, in [START_REF] Grégoire | Tracking mus and strict inconsistent cover[END_REF] the same authors have introduced a concept of inconsistent covers to circumvent the possible intractable number of MUSes, and presented a technique to compute them. Roughly, an inconsistent cover of an unsatisfiable SAT instance represents a set of MUSes that covers enough independent causes of inconsistency that would allow the instance to regain consistency if they were repaired. Although an inconsistent cover does not provide us with the set of all MUSes that may be present in a formula, it does however provide us with a series of minimal explanations of inconsistency that are sufficient to explain and potentially "fix" enough causes of inconsistency in order for the whole instance to regain consistency. These latter techniques are incomplete ones in the sense that they do not necessarily deliver all MUSes. However, in some application domains, it can be necessary to find the set of all MUSes, because diagnosing infeasibility is hard, if not impossible, without a complete view of its causes [START_REF] Liffiton | On finding all minimally unsatisfiable subformulas[END_REF]. Obviously enough, such techniques can only remain tractable provided that the number of MUSes remains itself tractable. Likewise, the number of MSSes and CoMSSes can be exponential in the worst case. It should be noted that many domains in Artificial Intelligence like belief revision (see e.g. [START_REF] Bessant | Iterated Syntax-Based Revision in a Nonmonotonic Setting[END_REF]) involve conceptual approaches to handle unsatisfiability that can require the complete sets of MUSes, MSSes, and CoMSSes to be computed in the worst case, even when additional epistemological ingredients like e.g. stratification are introduced in the logical framework.

In this paper, the focus is on complete techniques. We introduce a new complete technique to compute all MUSes, MSSes and CoMSSes of a SAT instance, provided obvious tractability limitations. It improves the currently most efficient complete technique, namely Liffiton and Sakallah'one [START_REF] Liffiton | On finding all minimally unsatisfiable subformulas[END_REF] (in short L&S), which in turn was shown more competitive than previous approaches by Bailey and Stuckey [START_REF] Bailey | Discovery of minimal unsatisfiable subsets of constraints using hitting set dualization[END_REF], and by de la Banda, Stuckey and Wazny [START_REF] De La Banda | Finding all minimal unsatisfiable subsets[END_REF], which were introduced in somewhat different contexts.

Our approach exhibits two main features. First, it is a hybridization of the L&S complete approach with a local search pretreatment. A local search technique is indeed used as an oracle to find potential CoMSSes of the SAT instance, which are themselves hitting sets of MUSes. We show that such a hybridization can yield exponential efficiency gains. Second, the efficiency of the approach relies on the crucial concept of critical clause, which appears to be a powerful ingredient of our technique to locate MUSes.

The rest of the paper is organized as follows. First, the reader is provided with the necessary background about SAT, MUSes and the dual concepts of MSSes and CoMSSes. Then, Liffiton and Sakallah's exhaustive approach is briefly presented. In Section 4, we show how this technique can be valuably hybridized with a local search pretreatment, making use of the critical clause concept. It is shown how this pretreatment can be theoretically valuable from a computational point of view. In Section 5, we compare this new approach with Liffiton and Sakallah's one on various benchmarks.

Background

In this section, the reader is provided with basic notions about SAT, MUSes, MSSes and CoMSSes.

Let L be a standard Boolean logical language built on a finite set of Boolean variables, noted a, b, c, etc. The ∧, ∨, ¬ and ⇒ symbols represent the standard conjunctive, disjunctive, negation and material implication connectives, respectively. Formulas and clauses will be noted using upper-case letters such as C. Sets of formulas will be represented using Greek letters like Γ or Σ. An interpretation is a truth assignment function that assigns values from {true, false} to every Boolean variable. A formula is satisfiable when there is at least one interpretation (called model) that satisfies it, i.e. that makes it become true. An interpretation will be noted by upper-case letters like I and will be represented by the set of literals that it satisfies. Actually, any formula in L can be represented (while preserving satisfiability) using a set (interpreted as a conjunction) of clauses, where a clause is a finite disjunction of literals, where a literal is a Boolean variable that is possibly negated. SAT is the NP-complete problem that consists in checking whether a set of Boolean clauses is satisfiable or not, i.e. whether there exists an interpretation that satisfies all clauses in the set or not.

When a SAT instance is unsatisfiable, it exhibits at least one Minimally Unsatisfiable Subformula, in short one MUS.

Definition 1 A MUS Γ of a SAT instance Σ is a set of clauses s.t. (1) Γ ⊆ Σ (2) Γ is unsatisfiable (3) ∀∆ ⊂ Γ, ∆ is satisfiable Example 2 Let Σ = {a, ¬c, ¬b ∨ ¬a, b, ¬b ∨ c}. Σ exhibits two MUSes, namely {a, b, ¬b ∨ ¬a} and {¬c, b, ¬b ∨ c}.
A dual concept is the notion of Maximal Satisfiable Subset (MSS) of a SAT instance.

Definition 3 A MSS Γ of a SAT instance Σ is a set of clauses s.t. (1) Γ ⊆ Σ (2) Γ is satisfiable (3) ∀∆ ⊆ (Σ \ Γ) s.t. ∆ = ∅, Γ ∪ ∆ is unsatisfiable.
The set-theoretical complement of a MSS w.r.t. a SAT instance is called a CoMSS.

Definition 4

The CoMSS of a MSS Γ of a SAT instance Σ is given by Σ \ Γ Example 5 Let us consider the formula Σ from the previous example. Σ exhibits five MSSes: {b}, {¬c, a}, {¬c, ¬b ∨ c}, {¬b ∨ c, ¬b ∨ ¬a}, {¬c, ¬b ∨ ¬a} As shown by several authors [START_REF] Liffiton | On finding all minimally unsatisfiable subformulas[END_REF], these concepts are correlated. Mainly, a CoMSS contains at least one clause from each MUS. Actually, a CoMSS is an irreducible hitting set of the set of MUSes. In a dual way, every MUS of a SAT instance is an irreducible hitting set of the CoMSSes. Accordingly, as emphasized by [START_REF] Liffiton | On finding all minimally unsatisfiable subformulas[END_REF], although MINIMAL-HITTING-SET is a NP-hard problem, irreducibility is a less strict requirement than minimal cardinality. Actually, a MUS can be generated in polynomial time from the set of CoMSSes.

Liffiton and Sakallah's Exhaustive Approach

Liffiton and Sakallah's approach [START_REF] Liffiton | On finding all minimally unsatisfiable subformulas[END_REF] to compute all MUSes (in short L&S) is based on the strong duality between MUSes and MSSes. To the best of our knowledge, it is currently the most efficient one. First it computes all MSSes before it extracts the corresponding set of MUSes. Here, the focus is on L&S first step since we shall improve it and adopt the second step as such.

L&S is integrated with a modern SAT solver and takes advantage of it. Roughly, every ith clause C i = x 1 ∨ ... ∨ x m of the SAT instance is augmented with a negated clause selector variable y i to yield C ′ i = x 1 ∨ ... ∨ x m ∨ ¬y i . While solving these new clauses, assigning y i to false has the effect of disabling or removing C i from the instance. Accordingly, a MSS can be obtained by finding a satisfying assignment with a minimal number of y i variables assigned false. The algorithm makes use of a sliding objective approach allowing for an incremental search. A bound on the number of y i that may be assigned to false is set. For each value of the bound, starting at 0 and incrementing by 1, an exhaustive search is performed for all satisfiable assignments to the augmented formula C ′ i , which will find all CoMSSes having their size equal to the bound. Whenever one solution is found, it is recorded, and a corresponding clause forcing out that solution (and any supersets of it) is inserted. This blocking clause is a disjunction of the y i variables for the clauses in that CoMSS.

Before beginning the search with the next bound, the algorithm checks that the new instance augmented with all the blocking clauses is still satisfiable without any bound on the y i variables. If there is no such satisfying assignment, the entire set of CoMSSes has been found and the algorithm terminates.

The second part of the algorithm computes the complete set of MUSes from the set of CoMSSes in a direct way. The approach that we shall introduce will include this second step as such.

Local Search and Critical Clauses

In this section, it is shown how the aforementioned exhaustive search algorithm can be improved in a dramatic manner by hybridizing it with an initial local search step, which provides valuable oracles for the subsequent exhaustive search process. We shall call the new approach HYCAM (HYbridization for Computing All Muses).

First, let us motivate our approach in an intuitive manner. Clearly, a (fast) initial local search run for satisfiability on the initial instance might encounter some actual MSSes. Whenever this phenomenon happens, it can prove valuable to record the corresponding CoMSSes in order to avoid computing them during the subsequent exhaustive search. Moreover, rather than checking whether we are faced with an actual MSS or not, it can prove useful to record the corresponding candidate CoMSS that will be checked later during the exhaustive search. Obviously enough, we must study which interpretations encountered during the local search process yield candidate MSSes and criteria must be defined in order to record a limited number of potentially candidate CoMSSes only. In this respect, a concept of critical clause will prove extremely valuable in the sense that it allows us to state necessary conditions for being a CoMSS that can be checked quickly. When all the remaining candidate CoMSSes are recorded, the incremental approach by Liffiton and Sakallah allows us to exploit this information in a very valuable and efficient way. Let us describe this in a more detailed manner.

A local search algorithm is thus run on the initial SAT instance. The goal is to record as many candidate CoMSSes as possible, based on the intuitive heuristics that local search often converges towards local minima, which could translate possibly good approximations of MSSes. A straightforward approach would consist in recording for each visited interpretation the set of unsatisfied clauses. Obviously enough, we do not need to record supersets of already recorded candidate CoMSSes since they cannot be actual CoMSSes as they are not minimal with respect to set-theoretic inclusion. More generally, we have adapted the technique proposed by Zhang in [START_REF] Zhang | On subsumption removal and on-the-fly cnf simplification[END_REF] to sets of clauses in order to record the currently smaller candidates CoMSSes among the already encountered series of sets of unsatisfied clauses. Now, crucial ingredients in our approach are the concepts of once-satisfied and critical clauses. Moreover, we have also exploited the following concept of critical clause, which has already proved valuable for locating MUSes and inconsistent covers using an incomplete technique based on local search [START_REF] Grégoire | Extracting MUS[END_REF] [START_REF] Grégoire | Tracking mus and strict inconsistent cover[END_REF].

Definition 6 A clause C is once-satisfied by an interpretation I iff exactly one literal of C is satisfied by I. A clause C that is falsified by the interpretation I is critical w.r.t. I iff the opposite of each literal of C belongs to at least one oncesatisfied clause by I.

Intuitively, a critical clause is thus a falsified clause that requires at least another

PROOF.

By definition, when a clause C f of Γ is not critical w.r.t. I, there exists at least one literal c ∈ C f whose truth-value can be inversed (i.e. flipped) without falsifying any other clause of Σ. Accordingly, Γ is not minimal and cannot be a CoMSS of Σ.

In practice, testing whether all falsified clauses are critical or not can be performed quickly and prevents many sets of clauses to be recorded as candidate CoMSSes.

Using these features, the local search run on the initial SAT instance, as described in Algorithm 1, yields a series of candidate CoMSSes. This information proves valuable and allows us to boost L&S complete search.

L&S is incremental in the sense that it computes CoMSSes of increasing sizes, progressively. After n iterations have been performed, all CoMSSes of cardinality lower or equal than n have been obtained. Accordingly, if we have recorded candidate CoMSSes containing n + 1 clauses, and if they are not supersets of already obtained CoMSSes, we are sure that they are actual CoMSSes. In this respect, we do not need to search them, and their corresponding blocking clauses can be inserted directly. Moreover, we do not need to perform the SAT test at the end of the n-th iteration, since we are then aware of the existence of larger CoMSSes.

It is also easy to show that the insertion of these blocking clauses can allow both NP-complete and CoNP-complete tests to be avoided. Let us illustrate this on an example.

Example 8 Let Σ be the following SAT instance. If the branching variables are chosen based on the lexical order, then a and b are assigned to true and c 5 is falsified. Thus, L&S tries to prove that this clause forms a CoMSS, which requires a NP-complete test (because it has to find a model of Σ\{¬a ∨ ¬b} ∪ {a, b}). On the contrary, when HYCAM is run, the blocking clause y 5 is added before the first iteration of the complete algorithm is performed, since y 5 is the clause selector variable of c 5 and since the local search has already delivered this CoMSS. In consequence, when a and b are assigned true, the DPLL-algorithm backtracks immediately as the {y 5 , ¬y 5 } unsatisfiable set has been obtained, without requiring any further NP-complete test.

Σ =                                                                                      c 0 : (d) c 1 : (b ∨ c) c 2 : (a ∨ b) c 3 : (a ∨ ¬c)
Similarly, the introduction of additional clause selector variables by HYCAM can reduce the number of CoNP-complete tests. For example, let us assume that e is the first branching variable, that e is assigned f alse and that the next variables are selected according to the lexical order. When a and b are assigned true, L&S tries to prove that {c 5 } is a CoMSS. Since ¬e is tautological consequence of Σ\{¬a ∨ ¬b} ∪ {a, b}, no model exists for Σ\{¬a ∨ ¬b} ∪ {a, b, ¬e}. Clearly, such a test is in CoNP. Thanks to the previously delivered candidate CoMSSes, HYCAM avoids this part of the search space to be explored. Indeed, since we know that {c 5 } is a CoMSS, when a and b are assigned true, no further CoNP tests are performed with respect to this partial assignment.

In fact, from a computational point of vue, the preliminary non-expensive local search eliminates nodes in the search tree, avoiding both NP and CoNP tests.

The HYCAM algorithm is depicted in Algorithm 2.

Experimental Evaluation

HYCAM has been implemented and compared to L&S from a practical point of view. For both algorithms, the complete search step is based on the use of MiniSat [START_REF] Eén | Minisat home[END_REF], which is currently one of the best modern SAT solvers. As a case study, we used Walksat [START_REF] Kautz | Walksat in the SAT'2004 competition[END_REF] for the local search pretreatment. The number of flips and tries of Walksat is related to the number of candidate CoMSS already found. For each try, a small number of flips is performed. If no new candidate is found during a try then a counter is incremented. When this counter exceeds a threshold (experimentally Algorithm 2: The HYCAM algorithm

Input: a CNF formula Σ Output: All MSSes of Σ begin cand ←-LS approximation(Σ) ; / * algorithm 1 * / Σ y ←-addSelectorClauses(Σ) ; k ←-0 ; M SS ←-∅ ; while SAT(Σ y ) do removeAllSetImplied({Σ\C|C ∈ M SS},cand) ; Σ y ←-addBlockingClausesOfSize(k,cand) ; M SS ←-M SS ∪ {Σ\C|C ∈ cand and |C| = k} ; M SS ←-M SS ∪ SAT with bound(k,Σ y ) ; k ←-k + 1
; return M SS ; end set to 30), we consider that no new candidate could be found by the local search. This way to end the local search pretreatment offers a good trade-off between the number of candidates found and the time spent. Besides, for all experiments, the time consumed by the local search step was less than 5% of the global time. All our experimental studies have been conducted on Intel Xeon 3GHz under Linux CentOS 4.1. (kernel 2.6.9) with a RAM memory size of 2Go. In the following, a time-out limit has been set to 3 CPU hours.

First, in Table 1, we report experimental results about the computation of MSSes on pigeon-hole and xor-chains benchmarks [START_REF] Satlib | [END_REF], which are globally unsatisfiable in the sense that removing any one of their clauses makes the instance become satisfiable. Obviously enough, such instances exhibit a number of CoMSSes equals to their number of clauses, and the size of any of their CoMSS is one. A significant time gap can be observed in favor of HYCAM. The efficiency gain ratio is even more significant when the size of the instance increased. For these instances, the local search run often succeeds in finding all CoMSSes, and the complete step often reduces to an unsatisfiability test. On the contrary, L&S explores many more nodes in the search space to deliver the CoMSSes.

In Table 2, experimental results on more difficult benchmarks from the annual SAT competition [START_REF] Satlib | [END_REF] are described. Their number of MSSes is often exponential, and computing them often remains intractable. Accordingly, we have limited the search to CoMSSes of restricted sizes, namely we have set a size limit to 5 clauses. As our experimental results illustrate, HYCAM outperforms L&S. For example, let us consider rand net40-30-10. This instance contains 5831 MSSes (with the size of their corresponding CoMSSes less than 5). L&S and HYCAM deliver this result in 1748 and 386 seconds, respectively. For the ca256 instance, HYCAM has extracted 9882 MSSes in less than 5 minutes whereas L&S did not manage to Computing all MSSes, CoMSSes and MUSes are highly intractable issues in the worst case. However, it can make sense to attempt to compute them for some reallife applications. In this paper, we have improved the currently most efficient exhaustive technique, namely Liffiton and Sakallah's method, in several ways. Our experimental results show dramatic efficiency gains for MSSes, CoMSSes and MUSes extracting. One interesting feature of the approach lies in its anytime character for computing MSSes. MSSes of increasing sizes are computed gradually. Accordingly, we can put a bound on the maximum size of the CoMSSes to be extracted, limiting the computing resources needed to extract them. To some extent, both L&S and HYCAM prove more adapted to extract complete sets of MSSes and CoMSSes than complete sets of MUSes. Indeed, the procedure involves computing MSSes (and thus CoMSSes) first. In this respect, we agree with Liffiton and Sakallah that an interesting path for future research concerns the study of how MUSes could be computed progressively from the growing set of extracted MSSes.

Many artificial intelligence research areas have studied various problems involving the manipulation of MUSes, MSSes and CoMSSes, like model-based diagnosis, belief revision, inconsistency handling in knowledge and belief bases, etc. These studies are often conducted from a conceptual point of view, or from a worst-case complexity point of view, only. We believe that the practical computational progresses as such as the ones obtained in this paper can prove valuable in handling these problems practically. In this respect, future research could concentrate on deriving specific algorithms for these AI issues, exploiting results like the ones described in this paper.

c 4 :

 4 (¬b ∨ ¬e) c 5 : (¬a ∨ ¬b) c 6 : (a ∨ e) c 7 : (¬a ∨ ¬e) c 8 : (b ∨ e) c 9 : (¬a ∨ b ∨ ¬c) c 10 : (¬a ∨ b ∨ ¬d) c 11 : (a ∨ ¬b ∨ c) c 12 : (a ∨ ¬b ∨ ¬d) Σ is an unsatisfiable SAT instance made of 13 clauses and making use of 5 variables. It exhibits 3 MUSes, which are illustrated in Figure 1, and admits 19 MSSes. Assume that both L&S and HYCAM are run on this instance. Its clauses are augmented by the ¬y i negated clause selector variables. Assume also that the local search performed by HYCAM provides 4 candidate CoMSSes: {c 5 }, {c 3 , c 2 }, {c 0 , c 1 , c 2 } and {c 3 , c 8 , c 10 }.

c 3 :Fig. 1 .

 31 Fig. 1. MUSes of Example 8

  Let Σ be a SAT instance and let I be an interpretation. Let Γ be a non-empty subset of Σ s.t. all clauses of Γ are all falsified by I. When at least one clause of Γ is not critical w.r.t. I, then Γ is not a CoMSS of Σ.

	Algorithm 1: Local Search approximation
	Input: a CNF formula Σ
	Output: Set of candidate CoMSS
	begin
	candidates ←-∅ ;
	#f ail ←-0 ;
	I ←-generate random interpretation() ;
	while (#f ail < PRESET #FAILURES AUTHORIZED) do
	newcandidates ←-FALSE ;
	for j = 1 to #FLIPS do
	Let ∆ be the set of falsified clauses by I ;
	if ∀C ∈ ∆, C is critical and ∆ is not implied in candidates then
	removeAllSetImplied(∆,candidates) ;
	candidates ←-∆ ∪ candidates ;
	newcandidates ←-TRUE ;
	flip(I) ;
	if not(newcandidates) then #f ail ←-#f ail + 1;
	return candidates ;
	end
	clause to be falsified in order to become satisfied, performing a flip. The follow-
	ing proposition shows how this concept allows us to eliminate wrong candidates
	CoMSSes.
	Proposition 7

Table 1

 1 L&S vs. HYCAM on globally unsatisfiable instances produce this result within 3 hours. Let us note that HYCAM also delivers CoMSSes made of 5 clauses after its computation is ended since we know that all sets of 5 falsified clauses recorded by the local search run and that are not supersets of the obtained smaller CoMSSs are actually also CoMSSes.In Table3, experimental results on hard instances to compute the complete set of MSSes and MUSes are reported. As explained above, both L&S and HYCAM approaches require all MSSes to be obtained before MUSes are computed. By allowing complete sets of MSSes to be delivered in a shorter time, HYCAM allows the complete set of MUSes to be computed for more instances and in a faster manner than L&S does. Obviously enough, when the number of MSSes or MUSes are exponential, computing and enumerating all of them remain intractable.For instance, L&S was unable to compute all MSSes of the php-012-011 instance within 3 hours CPU time, and could thus not discover its single MUS. HY-CAM extracted it in 2597 seconds. On all instances exhibiting unique or a nonexponential number of MUSes, HYCAM was clearly more efficient than L&S. For example, on the dlx2 aa instance, L&S and HYCAM discovered the 32 MUSes within 3.12 and 0.94 seconds, respectively. Let us note that the additional time spent to compute all MUSes from the set of MSSes is often very small unless of course the number of MUSes is exponential.

	Instance	(#var,#cla)	#MSSes	L&S (sec.)	HYCAM (sec.)
	hole6	(42,133)	133	0.040	0.051
	hole7	(56,204)	204	0.75	0.33
	hole8	(72,297)	297	33	1.60
	hole9	(90,415)	415	866	30
	hole10	(110,561)	561	7159	255
	x1 16	(46,122)	122	0.042	0.041
	x1 24	(70,186)	186	7.7	0.82
	x1 32	(94,250)	250	195	28
	x1 40	(118,314)	314	2722	486

6 Conclusions and Future Research

Table 3

 3 L&S vs. HYCAM on computing all MUSes

Acknowledgments

We thank Mark Liffiton for making his system available to us. This work has been supported in part by the Région Nord/Pas-de-Calais. A preliminary version of this paper was published in the proceedings of IJCAI'07.


Sylvain Pogodalla

Comments to Yao and Zou on: "Hybrid Categorial Type Logics and the Formal Treatment of Chinese"

Introduction

The paper "Hybrid Categorial Type Logics and the Formal Treatment of Chinese" presents the modeling of several linguistic phenomena occurring in Chinese. It relies on the formal framework of Hybrid Categorial Logics (HCL). This framework takes benefit from the ability to parametrize resource-sensitive logical systems with structural properties (such as associativity and commutativity), resulting in four systems: NL: non-associative and non-commutative Lambek calculus [START_REF] Lambek | On the calculus of syntactic types[END_REF]; NLP: non-associative and commutative Lambek calculus; L: associative and non-commutative Lambek calculus [START_REF] Lambek | The mathematics of sentence structure[END_REF]; LP: associative and commutative Lambek calculus, also known as (the intuitionistic fragment of) Linear Logic [START_REF] Girard | Linear logic[END_REF].

Each of the binary connectives of these systems can be indexed by a composition mode i to allow for the formulas 4 F where A is a set of atomic types:

F ::= A|F \ i F |F / i F |F • i F
The key point of hybrid systems is to offer a way to move from one deduction system to another with inclusion and interaction postulates [START_REF] Moortgat | Categorial type logics[END_REF] such as (1) or ( 2).

(1)

Γ[(∆ 1 , ∆ 2 ) • i ] ⇒ A Γ[(∆ 2 , ∆ 1 ) • i ] ⇒ A (2) Γ[((∆ 1 , ∆ 2 ) • i , ∆ 3 ) • j ] ⇒ A Γ[(∆ 1 , (∆ 2 , ∆ 3 ) • j ) • i ] ⇒ A
which allow for interaction between structures and formulas lying in different systems.

The authors emphasize that their modeling consider multimodality restricted to binary connectives, and that they favor not considering unary connectives as MultiModal Categorial Grammar (MMCG) [START_REF] Moortgat | Categorial type logics[END_REF] usually does. They show how this hybrid system allows them to give interesting accounts of several phenomena in Chinese. They focus in particular on various ways to express topicalization and conjunction of non constituents.

The modeling they propose brings to my attention different questions. Some of them, discussed in Section 2, relate to the formal system itself and to the motivations of using a multimodal system with only binary connectives rather than unary ones. They also relate to the design chosen for the interaction postulates. Another set of questions that Section 3 points out relate to the modeling capacity of the approach, and to the possible over-generation of the system.

2 On Formal Properties of Hybrid Categorial Grammars

Unary and Binary Operations

The paper puts a strong emphasis on the advantages of systems that use multimodality and structural postulates for binary connectives over those using unary connectives: "extensive use of structural modalities tends to result on very complex analyses" (in the abstract), "where [structural modalities] are used extensively, unduly complicated accounts tend to result" (end of the first section). However, we would like to have a better characterization of the nature of the involved complexity. Does it relate to the number of postulates? To the design of lexical categories and grammar engineering? To some computational issues?

It would be of a great value that the authors discuss the formal properties of the framework of hybrid logics with respect to MMCG. In particular, it is worth mentioning that there exists embedding theorems [START_REF] Hepple | Hybrid categorial logics[END_REF][START_REF] Kurtonina | Structural control[END_REF][START_REF] Moortgat | Categorial type logics[END_REF] between the different logics by mean of and ♦ unary connectives. These theorems use translations (•) from formulas of a strong logic L 1 into formulas of some weaker logic L 1 5 . They state that

L 1 A ⇒ B iff L 0 ♦ + R ♦ A ⇒ B
where L 0 ♦ is L 0 augmented with the unary connectives and R ♦ the relevant structural postulates. Would it be possible that any lexical item whose category requires some constraint relaxation and is thus expressed in a strong logic is modeled in the weaker logic L using the unary connectives. For instance, the topicalization of noun phrases of example (3.5) is modeled with s T /(s np) 6 using LP. Would it be possible to model this phenomenon using the unary connectives, possibly s T /( s/np) and the adequate translation of the structured antecedent?

Postulates

There also are some questions regarding the structural postulates and the interaction between the different composition modes. The paper should make them precise. The inclusion postulates are clearly stated. It is also clear which modes allow for permutation or associativity or combination of both. However, it should be explicitly stated that the associativity rule is bidirectional, as the examples show. That is both

Γ[((∆ 1 , ∆ 2 ) • , ∆ 3 ) • ] ⇒ A [A] Γ[(∆ 1 , (∆ 2 , ∆ 3 ) • ) • ] ⇒ A and Γ[(∆ 1 , (∆ 2 , ∆ 3 ) • ) • ] ⇒ A [A] Γ[((∆ 1 , ∆ 2 ) • , ∆ 3 ) • ] ⇒ A stand, what is usually written Γ[((∆ 1 , ∆ 2 ) • , ∆ 3 ) • ] ⇒ A [A] Γ[(∆ 1 , (∆ 2 , ∆ 3 ) • ) • ] ⇒ A
More importantly, the examples the paper provides don't make any use of the interaction postulate (rule (2.5) of the paper). Does it mean that it is not useful for the given examples but are useful for other phenomena? Does it mean the authors prefer not to consider this postulate? It would be interesting to comment on the scope of this postulate, both from a cross-linguistic perspective and from a Chinese language perspective. Also note that the interaction postulate slightly differs from the one given in [START_REF] Hepple | Hybrid categorial logics[END_REF]:

Γ[(∆ 1 , (∆ 2 , ∆ 3 ) • i ) • j ] ⇒ A Γ[((∆ 1 , ∆ 2 ) • i , ∆ 3 ) • j ] ⇒ A here vs. Γ[(∆ 1 , (∆ 2 , ∆ 3 ) • j ) • i ] ⇒ A Γ[((∆ 1 , ∆ 2 ) • i , ∆ 3 ) • j ] ⇒ A in [5]
3 Linguistic Modeling with Hybrid Categorial Grammars

The end of the paper is devoted to providing analysis of linguistic phenomena in Chinese. It focuses on topicalization and conjunction of non constituents.

Specificities of Chinese

While the modeling that are proposed here are convincing, it would help the reader to have a better exposition of the specificities of Chinese. Regarding topicalization, as far as I can see, only example (3.5) "topicalized complex head word of subject" really needs inclusion postulates. For instance in example 3.3, we can associate the type s T /((np/np)\s) to Zhe-jian instead of s T /((np/np) s) and still analyze Zhe-jian yifu buliao bu-cuo as the following derivation shows:

np ⇒ np np ⇒ np s ⇒ s [\L] (np, np\s) • ⇒ s [/L] ((np/np, np) • , np\s) • ⇒ s [A] (np/np, (np, np\s) • ) • ⇒ s [\R] (np, np\s) • ⇒ (np/np)\s s t ⇒ s T [/L] (s T /((np/np)\s), (np, np\s) • ) • ⇒ s T Zhe-jian yifu buliao bu-cuo
In order to require the commutative implication, the example should illustrate an extraction that is not peripheral, as in ( 4) where the cleft word relates to the object rather than to the subject, contrasting with (3). Note that with the type s T /((np/np) s) for zhe-ge would enforce the grammaticality of (4) 7 as the derivation [START_REF] Hepple | Hybrid categorial logics[END_REF] shows. 8 Such contrasts would help the reader without any knowledge of Chinese to grasp the insights of the modeling proposed here.

(

• ) • ⇒ s [A] ((np, (np\s)/np) • , np) • ⇒ s [/L] ((np, (np\s)/np) • , (np/np, np) • ) • ⇒ s [<] ((np, (np\s)/np) • , (np/np, np) ⊗ ) • ⇒ s [P] ((np, (np\s)/np) • , (np, np/np) ⊗ ) • ⇒ s [(2)] (((np, (np\s)/np) • , np) • , np/np) ⊗ ⇒ s [P] (np/np, ((np, (np\s)/np) • , np) • ) ⊗ ⇒ s [ R] ((np, (np\s)/np) • , np) • ⇒ (np/np) s s t ⇒ s T [/L] (s T /((np/np) s), ( (np, (np\s)
/np) • , np) • ) • ⇒ s T zhe-ge ban shang zhuxi xihuan xuesheng
A similar comment apply to example (3.4). The following derivation:

np\s ⇒ np\s np ⇒ np s ⇒ s [\L] (np, np\s) • ⇒ s [\L] (np, (np\s, (np\s)\(np\s)) • ) • ⇒ s [A] ((np, np\s) • , (np\s)\(np\s)) • ⇒ s [/R] (np, np\s) • ⇒ s/((np\s)\(np\s)) s T ⇒ s T [/L] (s T /(s/((np\s)\(np\s))), (np, np\s) • ) • ⇒ s T
shows that ( 6) is grammatical with this typing.9 [START_REF] Kurtonina | Structural control[END_REF] Finally, one could improve the understanding of the reader in providing glosses with a better alignment. It also seems that the examples include expressions that are provided with a single type, such as Zhe-jian : s T /((np/np) s) whereas it is itself a combination of words. It would be worth providing the very lexical types so that we can better understand the combination even if the derivations make use of partially evaluated expression. It amounts to ask: could (s T /((np/np) s))/n be a suitable type for zhe and n a suitable type for jian?

Conclusion

This paper provides an interesting use of hybrid systems to model Chinese. Defining more precisely the specificity the system should handle would improve the relevance of the approach. As such systems require fine-tuning to avoid under-and over-generation, automatic systems such as Grail [START_REF] Moot | Grail: an automated proof assistant for categorial grammar logics[END_REF] could be used.

3. 2

 2 Parsing, Proof-Search, and Over-generation My final comment has to do with choosing which sequent to prove. Indeed, as soon as several composition modes are available, one needs to specify the ones that are used in the structured antecedent of the sequent to prove. This choice is very important since choosing one composition mode or the other could prevent to use relevant structural postulates in the derivations and make them fail. So what are the strategies? Do we have to try all the possible composition modes? And what happens in case several structures would allow us to prove the sequent? It is indeed striking that examples (3.3) and (3.6) only make use of the L composition mode, while (3.4) and (3.5) mix L and LP composition modes, and (3.7) uses only LP composition mode. The latter is very surprising because it implies that the grammar will generate all the permutations: if (np , ((np\s)/np , ((np\s)\(np\s) , np) ⊗ ) ⊗ ) ⊗ ⇒ s Zhangsan chi le fan is derivable, so are Zhangsan chi fan le, Zhangsan le chi fan, Zhangsan le fan chi, Zhangsan fan chi le, Zhangsan fan le chi, chi Zhangsan le fan, chi Zhangsan fan le etc. I really wonder if all of them are acceptable.

  T /(s/np), ((np, (np\s)/np) • , np/np) • ) • ⇒ s T None of these examples makes use of an inclusion or interaction postulate. So it is really important to provide relevant examples focusing on the requirement Chinese makes on peripheral or non-peripheral extraction. The examples for conjunction don't make clear either what the specificities for Chinese are. If there are none, related works to conjunction for other languages should be mentioned.

	Ni-men-ban		we	shouxuan Yangfan
	You-PL-class		I	first-pick Yangfan
	s T /(s/((np\s)\(np\s)))	np	(np\s)
	And (7) illustrates a similar remark for example (3.6)
	(7) Jiu	wo	xihuan he	guizhou de
	Alcool	I	like drink	Guizhou NOMI
	s T /(s/np)	np	(np\s)/np	np/np
	(s

And the corresponding family of the structural connectives (•, •) • i that describes the structured antecedent in the sequent formulation of the calculus.

Different kinks of translations are considered in the above mentioned papers, but it seems to me it is the ones that are relevant here.

We follow the usual way to express the right implication of L with / and the implication of LP with rather than with indexed connectives.

To the best of my knowledge, this reading seems at least hard to get.

Note that the derivation uses interaction postulate (2) . It's possible that other examples don't. I let it to Chinese speakers!

I don't quite understand why the authors call this case "object modifier" since its rather a verb phrase modifier that is hypothesized in the derivation.


Pierre-Henri Conze

Tomás Crivelli

Philippe Robert

Luce Morin

⋆ Technicolor

DENSE MOTION ESTIMATION BETWEEN DISTANT FRAMES: COMBINATORIAL MULTI-STEP INTEGRATION AND STATISTICAL SELECTION

Keywords: motion estimation, statistical analysis, dense point matching, distant frames

Accurate estimation of dense point correspondences between two distant frames of a video sequence is a challenging task.

To address this problem, we present a combinatorial multistep integration procedure which allows one to obtain a large set of candidate motion fields between the two distant frames by considering multiple motion paths across the video sequence. Given this large candidate set, we propose to perform the optimal motion vector selection by combining a global optimization stage with a new statistical processing. Instead of considering a selection only based on intrinsic motion field quality and spatial regularization, the statistical processing exploits the spatial distribution of candidates and introduces an intra-candidate quality based on forward-backward consistency. Experiments evaluate the effectiveness of our method for distant motion estimation in the context of video editing.

INTRODUCTION

Despite rapid and significant progress since early formulations [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF][START_REF] Horn | Determining optical flow[END_REF], optical flow estimation still remains an open issue with crucial implications in computer vision. State-of-the-art methods [START_REF] Zach | A duality based approach for realtime TV-L1 optical flow[END_REF][START_REF] Steinbrucker | Large displacement optical flow computation without warping[END_REF][START_REF] Sun | Secrets of optical flow estimation and their principles[END_REF][START_REF] Sundaram | Dense point trajectories by GPU-accelerated large displacement optical flow[END_REF][START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF][START_REF] Xu | Motion detail preserving optical flow estimation[END_REF][START_REF] Tao | SimpleFlow: A non-iterative, sublinear optical flow algorithm[END_REF] have shown to be highly accurate for estimating dense motion fields between consecutive frames of a video sequence. However, they show limitations when applied to distant frames. The classical optical flow assumptions are not verified in this case, especially for difficult situations such as illumination changes, large motion, temporal occlusions, zooming, non-rigid deformations, low color contrast and transparency.

Direct matching between distant frames can be thus sensitive to ambiguous correspondences. An alternative consists in computing the long-range displacement through concatenation of elementary optical flow fields. This can be done by temporal integration, similarly to dense point tracking algorithms [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF]. However, even small errors in the input fields can lead to large drifts in the final motion field.

A first step towards accurate long-range dense correspondences is to combine numerous estimations from direct matching and temporal integration. Following a similar approach to that presented in [START_REF] Lempitsky | FusionFlow: Discrete-continuous optimization for optical flow estimation[END_REF], one can select for each pixel the optimal motion vector among a set of candidate motion fields based on intrinsic motion field quality and spatial regularization. A more sophisticated processing, described in [START_REF] Crivelli | From optical flow to dense long term correspondences[END_REF][START_REF] Crivelli | Multi-step flow fusion: Towards accurate and dense correspondences in long video shots[END_REF], consists in sequentially merging a set of concatenated multi-step motion fields at intermediate frames up to the target frame. However, in either case, the optimal motion vector selection strongly depends on the same optical flow assumptions that frequently fail between distant frames. This issue could be partially compensated by complexifying the matching criteria ad-infinitum, but an uncertainty component is always present. This argues in favor of a statistical processing which takes into account the random nature of these perturbations among a large set of dense motion fields.

In this direction, we propose two main contributions to address the dense matching problem between distant frames. Firstly, we present a combinatorial multi-step integration method which allows one to get a large set of motion fields between two distant frames by considering multiple motion paths across the sequence. Secondly, once this motion candidate construction is performed, we apply a new approach to select the optimal motion field based on statistics and spatial regularization. Results for motion estimation between distant frames in the context of video editing are presented.

MOTION CANDIDATE CONSTRUCTION

Let us consider a sequence of N +1 RGB images {I n } n∈[[0,...,N ]] and let I a and I b be two distant frames of this sequence with 0 ≤ a < b ≤ N . In this paper, we focus on the frame pair {I a , I b } and our goal is to accurately estimate a dense motion field between these two frames. In this section, we aim at generating multiple motion maps between I a and I b .

Direct matching with multiple optical flow estimators

A first approach for building multiple motion candidates consists in considering a direct motion estimation using different optical flow methods [START_REF] Lempitsky | FusionFlow: Discrete-continuous optimization for optical flow estimation[END_REF]. Even if the considered estimators may fail in some regions, the idea is to pool the strengths of each one. Furthermore, the same estimator can be used several times by modifying its parameter settings. In addition, we can derive from each version a parametric motion field for which motion is constrained by a global transformation [START_REF] Wills | A feature-based approach for determining dense long range correspondences[END_REF]. Direct matching is more ambiguous as the distance between I a and I b increases. Due to a large motion range, the motion of periodic color patterns or uniform areas may not be correctly estimated. This supports a motion field construction stage using concatenation of various optical flow fields.

Combinatorial multi-step integration

Let us describe the concept of motion path as an alternative to direct matching for obtaining a displacement map between I a and I b . A motion path is obtained through concatenation of elementary optical flow fields across the video sequence. It links each pixel x a of I a to a corresponding position in I b . Elementary optical flow fields can be computed between consecutive frames or with different frame steps [START_REF] Crivelli | From optical flow to dense long term correspondences[END_REF][START_REF] Crivelli | Multi-step flow fusion: Towards accurate and dense correspondences in long video shots[END_REF], i.e. with larger inter-frame distances. Let S n = {s 1 , s 2 , . . . , s Qn } ⊂ {1, . . . , N -n} be the set of Q n possible steps at instant n. This means that the set of optical flow fields {v n,n+s1 , v n,n+s2 , . . . , v n,n+sQ n } is available.

Our objective is to produce a large set of motion maps between I a and I b as to form a significative set of samples upon which a statistical processing would be meaningful and advantageous. Given this objective, we propose to initially generate all the possible step sequences (i.e. combinations of steps) in order to join I b from I a . Each step sequence defines a motion path. Let Γ a,b = {γ 0 , γ 1 , . . . , γ K-1 } be the set of K possible step sequences γ i between I a and I b . Γ a,b is computed by building a tree structure (Fig. 1) where each node corresponds to a motion field assigned to a given frame for a given step value (node value). In pratice, the construction of the tree is done recursively starting from I a : we create for each node as many children as the number of steps available at the current instant. A child node is not generated when I b has already been reached (therefore, the current node is considered as a leaf) or if I b is passed given the considered step. Finally, once the tree has been built, going from the root node to leaf nodes gives Γ a,b , the set of step sequences. For illustration, the tree in Fig. 1 indicates the four step sequences that can be generated going from I 0 to I 3 with steps 1, 2 and 3:

Γ 0,3 = {{1, 1, 1}, {1, 2}, {2, 1}, {3}}.
Once all the possible step sequences γ i ∀i ∈ [[0, . . . , K -1]] between I a and I b have been generated, the corresponding Fig. 2: Generation of multi-step motion paths. For each pixel x 0 of I 0 , this gives a set of candidate positions in I 3 . motion paths can be constructed through 1st-order Euler integration. Starting from each pixel x a of I a and for each step sequence, this integration performs the accumulation of optical flow fields following the steps which form the current step sequence. Thus, with steps 1, 2 and 3, Fig. 2 illustrates the construction of the four possible motion paths (one for each step sequence of Γ 0,3 ) between I 0 and I 3 . Let f i j = a + j k=0 s i k be the current frame number during the construction of motion path i from I a where j is the step index within the step sequence γ i . For each γ i ∈ Γ a,b and for each step s i j ∈ γ i , we start from x a in order to iteratively compute:

x f i j = x f i j-1 + v f i j-1 ,f i j (x f i j-1 ) (1)
Once all the steps s i j ∈ γ i have been run through, we obtain x i b , the corresponding position in I b of x a of I a obtained with step sequence γ i . By considering all the step sequences, we finally get a large set of candidate positions in I b and this for each pixel x a of I a . Note that the occlusion maps attached to input motion fields are used to possibly stop the motion path construction. Considering an intermediate point x f i j during the construction, a step can be added only if the closest pixel to x f i j is considered as un-occluded for this step. Otherwise, the motion path is removed. In the following, the large set of candidate positions in I b is defined as

T a,b (x a ) = {x i b } ∀i ∈ [[0, . . . , K x a -1]] where K x a is the cardinal of T a,b (x a ).
Up to now, we have considered an exhaustive generation of step sequences for clarity. However, for very distant frames and for a large set of steps, it is not possible to consider all possible step sequences (computational and memory issues). For instance, for a distance of 30 frames and with steps 1, 2, 5 and 10, the number of possible motion paths is 5877241. Therefore, the procedure described above is performed on a reasonable number of step sequences and not for all as previouly assumed. Firstly, we limit the number of elementary vectors composing the motion paths by providing a maximum number of concatenations N c . Indeed, the concatenation of numerous vectors may lead to an important drift. Secondly, we randomly select N s motion paths among the remaining motion paths (N s determined by storage capacity). This selection is guided by the fact that the candidate vectors should not be highly correlated. The frequency of appearance of a given step at a given frame must be uniform among all the possible steps arising from this frame in order to avoid a systematic bias towards the more populated branches of the tree.

MOTION VECTOR SELECTION ON LARGE SETS

Statistical processing for motion vector selection

Given T a,b (x a ) = {x i b } i∈[[0,...,Kx a -1]]
, the set of candidate positions in I b obtained for each pixel x a of I a from the motion field construction stage described above, the objective is to select the optimal candidate position x * b ∈ T a,b (x a ) by exploiting the statistical information on the point distribution and the quality of each candidate. The idea is to assume a Gaussian model for the distribution of T a,b (x a ) and try to find its central value, x * b . Using the maximum likelihood estimator (MLE) and imposing the selection among elements of T a,b (x a ), the choice of the optimal candidate position x * b is defined by:

x * b = arg min x i b Kx a -1 j=0 j =i x j b -x i b 2 2 (2)
The assumption of Gaussianity can be largely perturbed by outliers. Consequently, a robust estimation of the distribution central value is necessary:

x * b = arg min x i b med j =i x j b -x i b 2 2 (3)
Finally, each candidate position x i b receives a corresponding quality score Q(x i b ) computed using the inconsistency values Inc(x i b ). Inc(x i b ) corresponds to the Euclidean distance to the nearest reverse (resp. direct) candidate among the distribution if x i b is direct (resp. reverse). We aim at assigning high quality to candidates for which the corresponding motion field between I a and I b is consistent with a motion 3). Low PSNR for first pairs are due to the foreground object which degrades the estimation.

field between I b and I a . Quality scores Q(x i b ) are computed as follows: the maximum and minimum values of Inc(x i b ) among all candidates are mapped from 0 to a predefined integer Q max . Intermediate inconsistency values are mapped to the line defined by these two values and the result is rounded to the nearest integer:

Q(x i b ) ∈ [[0, . . . , Q max ]]. The higher Q(x i b ), the smaller Inc(x i b )
. We aim at promoting candidates in the neighborhood of high quality candidates. In practice, Q(x i b ) is used as a voting mechanism [START_REF] Yin | Circuits and Systems II: Analog and Digital Signal Processing[END_REF]: while computing the medians in Equation ( 3), each sample x j b is considered Q(x j b ) times to set the occurrence of elements

x j b -x i b 2 2
which enforces the forward-backward motion consistency. The statistical processing being applied for each pixel independently, we describe in what follows a global optimization method which includes regularization.

Global optimization for motion vector selection

We perform a global optimization stage that fuses for each pixel motion candidates into a single optimal motion field, following the approach of [START_REF] Lempitsky | FusionFlow: Discrete-continuous optimization for optical flow estimation[END_REF]. We introduce L = {l xa } as a labeling of pixels x a of I a where each label indicates one of the candidates of T a,b (x a ). Let d lx a a,b be the corresponding motion vectors of candidates of T a,b (x a ). We define the following energy and minimize it with fusion moves [START_REF] Lempitsky | FusionFlow: Discrete-continuous optimization for optical flow estimation[END_REF][START_REF] Lempitsky | Fusion moves for Markov random field optimization[END_REF].

E a,b (L) = xa ρ d (C(x a , d lx a a,b (x a )) + Inc(x a + d lx a a,b (x a ))) + <xa,y a > α x a ,y a .ρ r ( d lx a a,b (x a ) -d ly a a,b (y a ) 1 ) (4)
The data term involves the matching cost C(x a , d Table 2: Registration and PSNR assessment with: the combinatorial integration followed by the global optimization (GO); by the statistical processing combined to GO (SP+GO); the multi-step fusion (MSF) method [START_REF] Crivelli | Multi-step flow fusion: Towards accurate and dense correspondences in long video shots[END_REF]. PSNR scores are computed on the whole images of Newspaper (Fig. 4).

are described in [START_REF] Lempitsky | FusionFlow: Discrete-continuous optimization for optical flow estimation[END_REF]. Fusion moves algorithm fuses candidates pair by pair up to getting an optimal field d * a,b but its application to a large set is limited by the computational load.

Motion vector selection framework

We propose to combine statistical processing and the above global optimization stage to combine simultaneously information about the point distribution, a robust selection based on the intrinsic motion field quality and a spatial regularization. For each x a ∈ I a , the statistical processing is applied to the whole set T a,b (x a ). Then, we select the N opt best candidates of the distribution with the criterion of median minimization of (3). Finally, fusion moves algorithm fuses by pairs these N opt candidates up to obtaining the best one.

RESULTS

Our experiments focus on frame pairs taken from three sequences: MPI S1 [START_REF] Granados | MPI-S1[END_REF], Hope and Newspaper. For the selected pairs, the combinatorial multi-step integration has been performed taking input elementary flow fields estimated with a 2D version of the disparity estimator of [START_REF] Robert | Disparity-compensated view synthesis for s3d content correction[END_REF]. For all the experiments, the parameters are : N c = 7, N s = 100, Q max = 2, N opt = 3. Steps 1, 2, 3, 4, 5, 15 and 30 have been used.

After this motion field construction stage, we have compared three selection procedures: 1) the statistical processing (SP), 2) the global optimization (GO) method solved by fu-sion moves [START_REF] Lempitsky | Fusion moves for Markov random field optimization[END_REF], 3) the statistical processing combined with the global optimization (SP+GO). Firstly, the final fields of each method have been compared through registration and PSNR assessment. For a given pair {I a , I b }, the final fields are used to reconstruct I a from I b through motion compensation and color PSNR scores are computed between I a and the registred frame for non-occluded pixels. Tables 1 and2 show quantitative comparisons through PSNR computed for various distances between I a and I b respectively on the kiosk of MPI S1 and on whole images of Newspaper. An example of registration of the kiosk for a distance of 21 frames is provided Fig. 3. Results show that SP is better than GO for all pairs. The low diversity of candidates at the output of SP limits the effect of regularization and explains the slight improvement between SP and SP+GO. The example of Fig. 3 is interesting due to the temporary occlusion of the kiosk which is jumped by multi-step motion paths. For this complex situation, SP+GO is more adapted than GO. Secondly, in the context of video editing, we evaluate the accuracy of SP+GO and GO by motion compensating in I b logos manually inserted in I a . Fig. 4 presents results for Hope and Newspaper with a distance of 18 and 70 frames respectively. For both cases, SP+GO shows a clear improvement compared to GO.

The proposed combinatorial integration combined to SP+GO gives better performance compared to the multistep fusion (MSF) method [START_REF] Crivelli | Multi-step flow fusion: Towards accurate and dense correspondences in long video shots[END_REF] according to PSNR scores of Table 2. The MSF method itself has been shown in [START_REF] Crivelli | Multi-step flow fusion: Towards accurate and dense correspondences in long video shots[END_REF] to outperform state-of-the-art methods such as [START_REF] Zach | A duality based approach for realtime TV-L1 optical flow[END_REF][START_REF] Sundaram | Dense point trajectories by GPU-accelerated large displacement optical flow[END_REF][START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF].

CONCLUSION

We perform long-term dense matching by considering multiple multi-step motion paths along the sequence. Given the resulting large set of motion candidates, we apply a selection procedure where the global optimization stage is preceded by a new statistical processing which exploits the spatial distribution and the intrinsic quality of candidates. It leads to better results compared to state-of-the-art methods.

Fig. 1 :

 1 Fig. 1: Generation of step sequences. Going from the root node to leaf nodes of this tree structure gives Γ a,b , the set of possible step sequences from I a to I b .

Fig. 3 :

 3 Fig.3: Source frames of the MPI S1 sequence[START_REF] Granados | MPI-S1[END_REF] and reconstruction of the kiosk of I 46 from I 25 with: 1) the statistical processing (SP), 2) the global optimization (GO) method solved by fusion moves[START_REF] Lempitsky | Fusion moves for Markov random field optimization[END_REF], 3) both combined (SP+GO).

2. 3 .

 3 Joint forward and backward processing Motion estimation can be enhanced by considering both forward and backward motion fields. Similarly to the forward direction, the set of backward motion fields from each pixel x b of I b to I a can be computed by considering multi-step backward motion paths. These backward motion fields can be inversed into forward motion fields in order to enrich T a,b (x a ). Thus, backward motion vectors from pixels of I b are projected into I a . For each one, we identify the nearest pixel of the arrival position. The corresponding vector from I b to I a is reversed and started from the previously identified nearest pixel which gives a new candidate for T a,b (x a ). Candidates of T a,b (x a ) which have been obtained through this procedure are defined as reverse. Otherwise, we call them direct.

  Frame pairs {25,45} {25,46} {25,47} {25,48} {25,49}

Fig. 4 :

 4 Fig. 4: a-d) Logo insertion in I 5036 and propagation to I 5054 (Hope). e-h) Logo insertion in I 230 and propagation to I 160 (Newspaper). We compare the global optimization (GO) method with the statistical processing (SP) combined to GO (SP+GO).

Table 1 :

 1 Comparison through registration and PSNR assessment between: 1) the statistical processing (SP), 2) the global optimization (GO), 3) SP+GO. PSNR scores are computed on the kiosk of MPI S1 (Fig.

	{25,50}


B Romaniuk
email: barbara.romaniuk@greyc.ismra.fr

M Desvignes
email: michel.desvignes@greyc.ismra.fr

M Revenu

M J Deshayes

Linear and Non-Linear Model for Statistical Localization of Landmarks

This paper presents and compares 3 methods for the statistical localization of partially occulted landmarks. In many real applications, some information is visible in images and some parts are missing or occulted. These parts are estimated by 3 statistical approaches : a rigid registration, a linear method derived from PCA, which represents spatial relationships, and a non linear model based upon Kernel PCA. Applied to the cephalometric problem, the best method exhibits a mean error of 3.3 mm, which is about 3 times the intra-expert variability.

Introduction

The goal of orthodontic and orthognatic therapy is to improve the interrelationships among craniofacial tissues. A cephalogram is a two-dimensional X-Ray image of the sagital skull projection [START_REF] Macnamara | A method of Cephalometric Evaluation[END_REF] [START_REF] Cretot | L'image téléradiographique en céphalometrie[END_REF]. It is used to evaluate these relationships. Cephalometric landmarks are bony landmarks and are first located on the radiograph. Distances and angles among these landmarks are compared with normative values to diagnose a patient's deviation from ideal form and to evaluate craniofacial growth characteristics, skeletal and dental disharmonies. It is also used to evaluate results and stability of various treatment approaches. This task is challenging and has been the subject of previous research [START_REF] Davis | A Blackboard Architecture for Automatic Cephalogram Analysis[END_REF] [START_REF] Hutton | An Evaluation of Active Shape Models for the Automatic Identification of Cephalometric Landmarks[END_REF]. Our goal is the realization of a computer vision system to obtain an objective and reproducible cephalometric analysis. Indeed, large inter-expert and intra-expert variability has been noticed [START_REF] Cretot | L'image téléradiographique en céphalometrie[END_REF]. The main source of errors is the precise identification of landmarks. The two main causes are the subjectivity in the interpretation of the landmark definitions and the positional repeatability of human experts. Landmarks are difficult to distinguish on images and interpretation needs a long training time. In a computerized method, the formal descriptions of landmarks used by clinicians are not directly transposable: we then use a statistical approach to provide an initial estimation of landmark positions, using statistical models and training sets.

During the past decade, there has been a lot of work in shape based approaches, for segmentation, registration or identification tasks. In fact, any application where the geometric comparison of objects is required needs a shape analysis. The pioneers of the subject of shape analysis are Kendall [START_REF] Kendall | Shape manifolds, procrustean metrics and complex projectives spaces[END_REF] and Bookstein [START_REF] Bookstein | On the cephalometrics of skeletal change[END_REF]. In these works, shape is defined as the remaining information after alignment (rotation, translation, scale) between two objects. In image analysis, Pentland [START_REF] Sclaroff | Modal Matching for Correspondence and Recognition[END_REF] has defined modal analysis and a similar idea has been used by Cootes [START_REF] Cootes | Active Appearance Models[END_REF] in the Active Shape Model (ASM) and Active Appearance Model (AAM). They both involve a Principal Component Analysis (PCA) to build a statistical shape model. In this model, the mean object and the variation around this mean position are both represented. AAM was used for cephalometric purpose by Hutton [START_REF] Hutton | An Evaluation of Active Shape Models for the Automatic Identification of Cephalometric Landmarks[END_REF] without sufficient accuracy.

Other methods related to this problem use elastic registration to align an image with a model. The model can be an image [START_REF] Woods | Automated Image Registration: I. General Methods and Intrasubject, Intramodality Validation[END_REF], an atlas [START_REF] Subsol | A general scheme for automatically building 3D morphometric anatomical atlases: application to a skull atlas[END_REF] or a set of landmarks [START_REF] Bookstein | Principal Warps : Thin-Plate Splines and the Decomposition of Deformations[END_REF]. Elastic registration is a powerful tool based upon physical models such as solid or fluid deformations and includes complex and non linear model. Yet, the variability of the shape is not represented. Some works on Kernel PCA [START_REF] Mika | Kernel PCA and De-Noising in Feature Spaces[END_REF] are very close to our method. Briefly, Kernel PCA maps the input data in a Feature Space (F-Space) using a non linear mapping. PCA is performed in the F-Space. The mean shape is given by the eigenvectors corresponding to the largest eigenvalues. In a classification problem, classification is done in the F-Space. In a localization problem, mean shape in the Fspace must be back-projected in the input space. The choice of the mapping and the back projection are difficult problems and are still open issues.

In this paper, we present a comparison of 3 methods to localize landmarks and their application to the cephalometric problem. The first method is a simple affine registration, the second method is based upon a linear PCA and the third one is a non linear method close to PCA.

Methods

In the cephalometric problem, orthodontists have annotated cephalograms with 14 landmarks (cephalometric points) on a training set of radiographs. We also use an a-priori knowledge, which is a common knowledge in all cephalometric analysis : there is an unknown spatial relation between the cranial contour and the cephalometric points. The main problem in the cephalometric analysis is to discover this relation. Fortunately, the cranial contour can be automatically detected and extracted from the image [START_REF] Desvignes | Détection et modélisation du contour crânien sur des radiographies[END_REF], and then sampled (16 points). Our training set of points is composed of the 14 cephalometric points and the sampled version of the cranial contour. From this data-set, a mean shape model is computed. To retrieve landmarks on a new image, the cranial contour is detected and sampled, cephalograms are registered and the mean shape model is used to estimate the position.

Linear affine model

The problem is to compute a mean shape from a training set of points. First, all the sets of the training base have to be aligned. Procustes Analysis is a common tool to register two sets. It is a one to one mapping. To avoid this mapping, we have approximated the cranial contour by an ellipse, with the following parameters :

x g , y g : center of ellipse, q : angle between first principal axis and Ox, a,b : length of the principal axis. The coordinates of cephalometric points are expressed in the coordinate space defined by the center x g , y g and the vectors a and b along the principal axis of the ellipse.

Let X i =(x i ,y i ) i Î1..n, be the points of the cranial contour. We can write :

å = = ÷ ÷ ø ö ç ç è ae = n i i g g g X n y x X 1 1 , å = - - = n i q g i p g i pq y y x x n 1 ) ( ) ( 1 m , 20 02 11 2 ) * 2 tan( m m m q - = .
Intersections between the cranial contour and the principal axis defined by the angles q and q+2p define the the unit vectors a and b.

Let C i j =(cx i j ,cy i j ), i Î1..m, jÎ1..P, be the i th cephalometric point of j th image. The mean shape C i of the i th cephalometric points is then defined by :

÷ ÷ ø ö ç ç è ae - - ÷ ÷ ø ö ç ç è ae ÷ ÷ ø ö ç ç è ae ÷ ÷ ø ö ç ç è ae - = - = å G j i G j i P j i y y x x b a P C * 0 0 * ) cos( ) sin( ) sin( ) cos( 1 ˆ1 1 q q q q .
On an unseen image, the cranial contour is detected and is fitted with an ellipse and the 5 parameters x g , y g, q , a and b are computed. The estimated landmarks are then :

( ) ÷ ÷ ø ö ç ç è ae + ÷ ÷ ø ö ç ç è ae ÷ ÷ ø ö ç ç è ae ÷ ÷ ø ö ç ç è ae - = ÷ ÷ ø ö ç ç è ae G G i i i y x C b a y x * 0 0 * ) cos( ) sin( ) sin( ) cos( q q q q .

Linear PCA model

In the previous method, spatial relations between cephalometric points are not examined although they seem to be quite important for the expert. The linear PCA method defined here is an elegant way to take into account spatial relations between landmarks and between landmarks and contour, and can also estimate the unknown part (cephalometric points) of the partially visible or occulted model (cranial contour). Let X i =(x 1i , y 1i , x 2i , y 2i , .. x ni , y ni ) ÎR 2n be the locations of the n cephalometric points on the i th cephalogram, C i be the locations of the m points of the sampled cranial contour on the i th cephalogram, and T i =(X i ,C i ) the concatenation of X i and C i . To compute a model with this training set, the first step is to align all these samples. This is realized with an iterative version of the Procuste analysis. Using PCA, we can write b T T i f + »

where :

å = = p i i T p T 1 1
is the mean shape of the pattern,

f = (f 1 | f 2 | f 3 |…| f t )
å = - - - = p i T i i T T T T p S 1 ) )( ( 1 1 , b is a vector of dimension t : ) ( T T b i T - = f
. The dimension t of the vector b is the number of eigenvectors with the largest eigenvalues. In classical use of PCA, t is chosen by

å å + = = ³ m n i i t i i 1 1 95 . 0 l l
, i.e. only eigenvectors that explain sufficiently the standard deviation are kept. The vector b of dimension t is a good approximation for the original data set and any set of n+m points can be represented or retrieved with the t (t<n+m) values of the vector b. Then, PCA can be seen as a denoising method. Now, on a new image, the cranial contour is extracted. To compare this contour with the model, we need to align it with the model. The transformation to align the n points (x i , y i ) with the model (x' i ,y' i ) is defined by a least square process :

÷ ÷ ø ö ç ç è ae + ÷ ÷ ø ö ç ç è ae ÷ ÷ ø ö ç ç è ae - = ÷ ÷ ø ö ç ç è ae y x t t y x y x T a b b a , 2 2 2 ' ' ' . 1 1 1 1 x x x y n x n y y n x x n i i i i i i = + + = å å å å a , 1051-4651/02 $17.00 (c) 2002 IEEE ' ' 2 2 ' ' 1 1 1 1 1 1 x x y n y x n y n x n x y n y x n i i i i i i i i i i å å å å å å - = + - = b .
Once the registration of a new form onto the model is done, the original idea of the paper is the following one : if some information is present, if some points (saying t=n points) are known, the remaining unknown m points can be determined using PCA, under the hypothesis that the first n eigenvalues of the covariance matrix explain the training dataset. Let C=(x,y) be the points of the cranial contour and X be the cephalometric points. Without any approximations, we can write :

÷ ÷ ÷ ÷ ÷ ÷ ÷ ø ö ç ç ç ç ç ç ç è ae ÷ ÷ ÷ ÷ ÷ ø ö ç ç ç ç ç è ae + ÷ ÷ ÷ ÷ ÷ ÷ ÷ ÷ ø ö ç ç ç ç ç ç ç ç è ae = ÷ ÷ ÷ ÷ ÷ ÷ ÷ ø ö ç ç ç ç ç ç ç è ae + + + + + + m n n m n m n m n m n m n m n m n b b b b b X X C C C X X C C C 1 2 1 , 2 , 1 , 1 , 2 , 1 2 , 1 1 , 1 1 2 1 1 2 1 .. .. M M M M M M M M M f f f f f f f
This is a linear system with n+m equations and n+2m unknown (X 1 ,X m ,b 1 ,..,b m+n ) that can not be resolved. Since PCA represents the dataset with less values, we can write, using t=n :

1 , 2 , 1 , 1 , 2 , 1 2 , 1 1 , 1 1 2 1 1 2 1 f f f f f f f M M M M M M M M M .
Now the unknown vector (b 1 ,..,b n ) is given by the n first equations of the system. Notice that if we choose t<n, the system is now overdetermined (t<n), a least square method is used to resolve the system. The m last equations are a linear system with m equations and m unknowns X 1 ..X m :

÷ ÷ ÷ ÷ ÷ ø ö ç ç ç ç ç è ae ÷ ÷ ÷ ÷ ÷ ø ö ç ç ç ç ç è ae - ÷ ÷ ÷ ÷ ø ö ç ç ç ç è ae ÷ ÷ ÷ ÷ ÷ ø ö ç ç ç ç ç è ae = ÷ ÷ ÷ ÷ ø ö ç ç ç ç è ae - n n n n n n n n n C C C C C C b b b M M L M O M M L L M 2 1 2 1 1 , 2 , 1 , , 2 2 , 2 1 , 2 , 1 2 , 1 1 , 1 2 1 f f f f f f f f f
+ + + + + + n n m n m n m n n n n n n m b b b X X M L M O M M L M 2 1 , 2 , 1 , , 1 2 , 1 , 1 1 f f f f f f .
In this framework, a linear approximation of spatial relations between cranial contour and cephalometric points are explicitly determined from the eigenvectors of the covariance matrix. In the same way, relations between landmarks are approximated.

Non Linear model

In this model, we project the location of the initial data in a new space, following the idea of Kernel methods .

Let P i be the points of the cranial contour. The origin of the sampling is given by the point P 0 with the higher curvature, i.e. a point near the intersection of the curve and of the nasal bone. Let the set of vectors R be defined by each pair of points P i and P j of the sample cranial contour.

(

)

( ) { } þ ý ü î í ì = = $ < Î " = ® k i i j i i i P P v and P O i k j n k j v O R r r / , , ... 1 , ,
Let M (M x ,M y ,1) a point. The coordinates a i of M in R are defined by the scalar product (or projection) of M and each vectors of R, i.e.

i v M i O i r ® = a .
Then, we can write :

() ( ) ( ) ÷ ÷ ÷ ø ö ç ç ç è ae ÷ ø ö ç è ae = ÷ ÷ ÷ ø ö ç ç ç è ae + - = = 1 1 A y x y x iy iy ix ix iy ix i M M M M v O v O v v a a .
The coordinates a i are the projection (mapping) of the point M on each vector built using the cranial contour. Our model represents the position and the variability of these new values. The mean E(a) and the standard deviation s(a) of the location of each landmark are computed and stored. As the number of coefficients a i is large n(n+1)/2, and because all the coefficients are not relevant for each landmark in the image, only the coefficients which are really useful (i.e. with small standard deviation) are stored.

To determine the unknown position of the landmark X in a new cephalogram, we can write:

[] X A E × = a where
A is computed by detecting the cranial contour and the vector a on the new cephalogram. X is the unknown cephalometric point and is computed with a weight least square inversion matrix, using the constraint X z =1.

Results

We have now 227 cephalograms annotated by an orthodontist. All the cephalograms have been processed by the 3 methods. The methods have failed on one cephalogram because the cranial contour was not detected properly. To test and compare the accuracy of these methods, this training set is divided in a first set of 60 cephalograms, randomly chosen, used to compute the model and a second set used to test the methods. We have measured the mean error between the real position and the estimated position of landmarks, on the x-axis and on the y-axis. Results are summarized in the table 1. For the cephalometric use, these values must be compared with intra and inter experts variability. For practical reasons, only the intra-expert variability have been estimated to 1.1mm on a small set (15) of cephalograms. In literature, inter-expert is often larger (5 times). Then, results of the statistical process are quite accurate, but not enough for clinical usage.

Conclusion

In this paper, we have addressed the problem of locating occulted parts of objects described by some landmarks. When only some of the landmarks are visible or extracted from an image, the remaining unknown landmarks must be estimated. A statistical method solves easily this problem. 3 methods have been described and compared.

Rigid registration or Procustes Analysis aligns (translation, rotation, scale) two or more sets of points, which allows to compute a mean model. This mean model is simply registered on a new image using the extracted (visible) landmarks.

Classical shape analysis uses PCA to determine a mean value and a set of allowed deformation modes and parameters. These modes are the eigenvectors of the covariance matrix, which described the spatial relationships between landmarks. The deformation parameters are computed from the extracted landmarks, and unknown landmarks are determined from these values. Spatial relationships and the shape of the actual image are then taken into account by this method.

At last, a non linear method maps the initial data in a large feature space using a projection. For each landmark, only the few useful axes of the space are used to estimate the unknown landmarks by a least square inversion.

Applied to the cephalometric problem, the set of landmarks is composed of cephalometric points (unknown) and of the cranial contour, which is automatically extracted. Results are quite satisfactory, but accuracy is not sufficient for clinical use.

Improvement of this accuracy is possible, and we are working on a non linear-method which is a mixture of the second and the third method : it will introduce spatial relationships in non-linear model and optimize the non linear mapping for the cephalometric problem.

2

 2

Figure 1 :

 1 Figure 1 : Annoted Points (black) and estimation (white) ãCranexplo

Table 1 :

 1

		Mean Error
	Error (mm) x	y
	Method 1	3.4	3.5
	Method 2	3.1	2.6
	Method 3	2.7	2.0

Acknowledgements

This work is achieved with Telecrane Innovation SA.