Guilhem Castagnos
email: guilhem.castagnos@math.u-bordeaux1.fr

Antoine Joux
email: antoine.joux@m4x.org

Fabien Laguillaumie

Phong Q Nguyen

Factoring pq 2 with Quadratic Forms: Nice Cryptanalyses

Keywords: Public-key Cryptanalysis, Factorisation, Binary Quadratic Forms, Homogeneous Coppersmith's Root Finding, Lattices

We present a new algorithm based on binary quadratic forms to factor integers of the form N = pq 2 . Its heuristic running time is exponential in the general case, but becomes polynomial when special (arithmetic) hints are available, which is exactly the case for the so-called NICE family of public-key cryptosystems based on quadratic fields introduced in the late 90s. Such cryptosystems come in two flavours, depending on whether the quadratic field is imaginary or real. Our factoring algorithm yields a general key-recovery polynomial-time attack on NICE, which works for both versions: Castagnos and Laguillaumie recently obtained a total break of imaginary-NICE, but their attack could not apply to real-NICE. Our algorithm is rather different from classical factoring algorithms: it combines Lagrange's reduction of quadratic forms with a provable variant of Coppersmith's lattice-based root finding algorithm for homogeneous polynomials. It is very efficient given either of the following arithmetic hints: the public key of imaginary-NICE, which provides an alternative to the CL attack; or the knowledge that the regulator of the quadratic field Q( √ p) is unusually small, just like in real-NICE.

Introduction

Many public-key cryptosystems require the hardness of factoring large integers of the special form N = pq 2 , such as Okamoto's Esign [START_REF] Okamoto | A fast signature scheme based on congruential polynomial operations[END_REF], Okamoto and Uchiyama's encryption [START_REF] Okamoto | A New Public-Key Cryptosystem as Secure as Factoring[END_REF], Takagi's fast RSA variants [START_REF] Takagi | Fast RSA-Type Cryptosystem Modulo p k q[END_REF], and the large family (surveyed in [START_REF] Buchmann | Number Field Cryptography[END_REF]) of cryptosystems based on quadratic fields, which was initiated by Buchmann and Williams' key exchange [START_REF] Buchmann | A Key-Exchange System based on Imaginary Quadratic Fields[END_REF], and which includes NICE 1 cryptosystems [HPT99,PT99,PT00,JSW08] (whose main feature is a quadratic decryption). These moduli are popular because they can lead to special functionalities (like homomorphic encryption) or improved efficiency (compared to RSA). And no significant weakness has been found compared to standard RSA moduli of the form N = pq: to the best of our knowledge, the only results on pq 2 factorisation are [START_REF] Peralta | Faster Factoring of Integers of a Special Form[END_REF][START_REF] Peralta | Elliptic curve factorization using a "partially oblivious" function. Cryptography and computational number theory[END_REF][START_REF] Boneh | Factoring N = p r q for large r[END_REF]. More precisely, [START_REF] Peralta | Faster Factoring of Integers of a Special Form[END_REF][START_REF] Peralta | Elliptic curve factorization using a "partially oblivious" function. Cryptography and computational number theory[END_REF] obtained a linear speed-up of Lenstra's ECM, and [BDH99, Sect. 6] can factor in time Õ(N 1/9 ) when p and q are balanced. Furthermore, computing the "squarefree part" of an integer (that is, given N ∈ N as input, compute (r, s) ∈ N 2 such that N = r 2 s with s squarefree) is a classical problem in algorithmic number theory (cf. [START_REF] Adleman | Open problems in number theoretic complexity[END_REF]), because it is polynomial-time equivalent to determining the ring of integers of a number field [START_REF] Chistov | The complexity of constructing the ring of integers of a global field[END_REF].

However, some of these cryptosystems actually provide additional information (other than N ) in the public key, which may render factorisation easy. For instance, Howgrave-Graham [START_REF] Howgrave-Graham | Approximate Integer Common Divisors Proc[END_REF] showed that the public key of [START_REF] Okamoto | Fast public-key cryptosystem using congruent polynomial equations[END_REF] disclosed the secret factorisation in polynomial time, using the gcd extension of Coppersmith's root finding method [START_REF] Coppersmith | Small Solutions to Polynomial Equations, and Low Exponent RSA Vulnerabilities[END_REF]. Very recently, Castagnos and Laguillaumie [START_REF] Castagnos | On the Security of Cryptosystems with Quadratic Decryption: The Nicest Cryptanalysis[END_REF] showed that the public key in the imaginary version [HPT99, PT99, PT00] of NICE allowed to retrieve the secret factorisation in polynomial time. And this additional information in the public key was crucial to make the complexity of decryption quadratic in imaginary-NICE, which was the main claimed benefit of NICE. But surprisingly, the attack of [START_REF] Castagnos | On the Security of Cryptosystems with Quadratic Decryption: The Nicest Cryptanalysis[END_REF] does not work against REAL-NICE [JSW08], which is the version of NICE with real (rather than imaginary) quadratic fields, and which also offers quadratic decryption. In particular, the public key of REAL-NICE only consists of N = pq 2 , but the prime p has special arithmetic properties. Our Results. We present a new algorithm to factor integers of the form N = pq 2 , based on binary quadratic forms (or equivalently, ideals of orders of quadratic number fields). In the worst case, its heuristic running time is exponential, namely Õ(p 1/2 ). But in the presence of special hints, it becomes heuristically polynomial. These hints are different from the usual ones of lattice-based factoring methods [START_REF] Coppersmith | Small Solutions to Polynomial Equations, and Low Exponent RSA Vulnerabilities[END_REF][START_REF] Boneh | Factoring N = p r q for large r[END_REF][START_REF] Howgrave-Graham | Approximate Integer Common Divisors Proc[END_REF] where they are a fraction of the bits of the secret prime factors. Instead, our hints are arithmetic, and correspond exactly to the situation of NICE, including both the imaginary [HPT99, PT99, PT00] and real versions [START_REF] Jacobson | An Adaptation of the NICE Cryptosystem to Real Quadratic Orders[END_REF]. This gives rise to the first general key-recovery polynomial-time attack on NICE, using only the public key.

More precisely, our arithmetic hints can be either of the following two:

i. The hint is an ideal equivalent to a secret ideal of norm q 2 in an imaginary quadratic field of discriminant -pq 2 : in NICE, such an ideal is disclosed by the public key. This gives an alternative attack of NICE, different from [START_REF] Castagnos | On the Security of Cryptosystems with Quadratic Decryption: The Nicest Cryptanalysis[END_REF]. ii. The hint is the knowledge that the regulator of the quadratic field Q( √ p) is unusually small, just like in REAL-NICE. Roughly speaking, the regulator is a

1 for New Ideal Coset Encryption real number which determines how "dense" the units of the ring of integers of the number field Q( √ p) are. This number is known to lie in the large interval log 1 2 ( √ p -4 + √ p) , 1 2 p 1 2 log p + 1 . But for infinitely many p (including square-free numbers of the form p = k 2 + r, where p > 5, r|4k and -k < r ≤ k, see [START_REF] Degert | Uber die Bestimmung der Grundeinheit gewisser reell-quadratischer Zhalkörper[END_REF]), the regulator is at most polynomial in log p. For these unusually small regulators, our algorithm heuristically runs in time polynomial in the bit-length of N = pq 2 , which gives the first total break of REAL-NICE [START_REF] Jacobson | An Adaptation of the NICE Cryptosystem to Real Quadratic Orders[END_REF]. We stress that although such p's are easy to construct, their density is believed to be arbitrary small. Interestingly, our algorithm is rather different from classical factoring algorithms. It is a combination of Lagrange's reduction of quadratic forms with a provable variant of Coppersmith's lattice-based root finding algorithm [START_REF] Coppersmith | Small Solutions to Polynomial Equations, and Low Exponent RSA Vulnerabilities[END_REF] for homogeneous polynomials. In a nutshell, our factoring method first looks for a reduced binary quadratic form f (x, y) = ax 2 + bxy + cy 2 representing properly q 2 with small coefficients, i.e. there exist small coprime integers x 0 and y 0 such that q 2 = f (x 0 , y 0 ). In case i., such a quadratic form is already given. In case ii., such a quadratic form is found by a walk along the principal cycle of the class group of discriminant pq 2 , using Lagrange's reduction of (indefinite) quadratic forms. Finally, the algorithm finds such small coprime integers x 0 and y 0 such that q 2 = f (x 0 , y 0 ), by using the fact that gcd(f (x 0 , y 0 ), pq 2 ) is large. This discloses q 2 and therefore the factorisation of N . In both cases, the search for x 0 and y 0 is done with a new rigorous homogeneous bivariate variant of Coppersmith's method, which might be of independent interest: by the way, it was pointed out to us that Bernstein [START_REF] Bernstein | List decoding for binary Goppa codes[END_REF] independently used a similar method in the different context of Goppa codes decoding.

Our algorithm requires "natural" bounds on the roots of reduced quadratic forms of a special shape. We are unable to prove rigorously all these bounds, which makes our algorithm heuristic (like many factoring algorithms). But we have performed many experiments supporting such bounds, and the algorithm works very well in practice. Factorisation and Quadratic Forms. Our algorithm is based on quadratic forms, which share a long history with factoring (see [START_REF] Crandall | Prime Numbers: A Computational Perspective[END_REF]). Fermat's factoring method represents N in two intrinsically different ways by the quadratic form x 2 + y 2 . It has been improved by Shanks with SQUFOF, whose complexity is Õ(N 1/4 ) (see [START_REF] Gower | Square form factorization[END_REF] for a detailed analysis). Like ours, this method works with the infrastructure of a class group of positive discriminant, but is different in spirit since it searches for an ambiguous form (after having found a square form), and does not focus on discriminants of a special shape. Schoof's factoring algorithms [START_REF] Schoof | Quadratic fields and factorization[END_REF] are also essentially looking for ambiguous forms. One is based on computation in class groups of complex quadratic orders and the other is close to SQUFOF since it works with real quadratic orders by computing a good approximation of the regulator to find an ambiguous form. Like SQUFOF, this algorithm does not takes advantage of working in a non-maximal order and is rather different from our algorithm. Both algorithms of [START_REF] Schoof | Quadratic fields and factorization[END_REF] runs in Õ(N 1/5 ) under the generalised Riemann hypothesis. McKee's method [START_REF] Mckee | Speeding Fermat's factoring method[END_REF] is a speedup of Fermat's algorithm (and was presented as an alternative to SQUFOF) with a heuristic complexity of Õ(N 1/4 ) instead of Õ(N 1/2 ).

SQUFOF and other exponential methods are often used to factor small numbers (say 50 to 100 bits), for instance in the post-sieving phase of the Number Field Sieve algorithm. Some interesting experimental comparisons can be found in [START_REF] Milan | Factoring Small Integers: An Experimental Comparison[END_REF]. Note that the currently fastest rigorous deterministic algorithm actually has exponential complexity: it is based on a polynomial evaluation method (for a polynomial of the form x(x -1) • • • (x -B + 1) for some bound B) and its best variant is described in [START_REF] Bostan | Linear Recurrences with Polynomial Coefficients and Application to Integer Factorization and Cartier-Manin Operator[END_REF]. Finally, all sieve factoring algorithms are somewhat related to quadratic forms, since their goal is to find random pairs (x, y) of integers such that x 2 ≡ y 2 mod N . However, these algorithms factor generic numbers and have a subexponential complexity. Road Map. The rest of the paper is organised as follows. The first section recalls facts on quadratic fields and quadratic forms, and present our heuristic supported by experiments. The next section describes the homogeneous Coppersmith method and the following exhibits our main result: the factoring algorithm. The last section consists of the two cryptanalyses of cryptosystems based on real quadratic fields (REAL-NICE) and on imaginary quadratic fields (NICE).

2 Background on Quadratic Fields and Quadratic Forms

Quadratic Fields

Let D = 0, 1 be a squarefree integer and consider the quadratic number field K = Q( √ D). If D < 0 (resp. D > 0), K is called an imaginary (resp. a real ) quadratic field. The fundamental discriminant ∆ K of K is defined as ∆ K = D if D ≡ 1 (mod 4) and ∆ K = 4D otherwise. An order O in K is a subset of K such that O is a subring of K containing 1 and O is a free Z-module of rank 2. The ring O ∆ K of algebraic integers in K is the maximal order of K. It can be written as Z + ω K Z, where

ω K = 1 2 (∆ K + √ ∆ K ). If we set f = [O ∆ K : O] the finite index of any order O in O ∆ K , then O = Z + f ω K Z. The integer f is called the conductor of O. The discriminant of O is then ∆ f = f 2 ∆ K . Now, let O ∆ be an order of discriminant ∆ and a be a nonzero ideal of O ∆ , its norm is N(a) = |O ∆ /a|. A fractional ideal is a subset a ⊂ K such that da is an ideal of O ∆ for d ∈ N. A fractional ideal a is said to be invertible if there exists an another fractional ideal b such that ab = O ∆ . The ideal class group of O ∆ is C(O ∆ ) = I(O ∆ )/P (O ∆ ), where I(O ∆ ) is the group of invertible fractional ideals of O ∆ and P (O ∆ ) the subgroup consisting of principal ideals. Its cardinality is the class number of O ∆ denoted by h(O ∆ ). A nonzero ideal a of O ∆ , a is said to be prime to f if a + f O ∆ = O ∆ . We denote by I(O ∆ , f ) the subgroup of I(O ∆ ) of ideals prime to f . The group O ⋆ ∆ of units in O ∆ is equal to {±1} for all ∆ < 0, except when ∆ is equal to -3 and -4 (O ⋆
-3 and O ⋆ -4 are respectively the group of sixth and fourth roots of unity). When ∆ > 0, then

O ⋆ ∆ = -1, ε ∆ where ε ∆ > 0 is called the fundamental unit. The real number R ∆ = log(ε ∆ ) is the regulator of O ∆ .
The following important bounds on the regulator of a real quadratic field can be found in [START_REF] Jacobson | An investigation of bounds for the regulator of quadratic fields[END_REF]:

log 1 2 ( √ ∆ -4 + √ ∆) ≤ R ∆ < 1 2 ∆ 1 2 log ∆ + 1 . (1)
The lower bound is reached infinitely often, for instance with ∆ = x 2 + 4 with 2 ∤ x. Finally, this last proposition is the heart of both NICE and REAL-NICE.

Proposition 1 ([Cox99, Proposition 7.20] [Wei04, Theorem 2.16]). Let O ∆ f be an order of conductor f in a quadratic field K. i. If A is an O ∆ K -ideal prime to f , then A ∩ O ∆ f is an O ∆ f -ideal prime to f of the same norm. ii. If a is an O ∆ f -ideal prime to f , then aO ∆ K is an O ∆ K -ideal prime to f of the same norm. iii. The map ϕ f : I(O ∆ f , f ) → I(O ∆ K , f ), a → aO ∆ K is an isomorphism. The map ϕ f from Proposition 1 induces a surjection φf : C(O ∆ f ) ։ C(O ∆ K )
which can be efficiently computed (see [START_REF] Paulus | A New Public-Key Cryptosystem over a Quadratic Order with Quadratic Decryption Time[END_REF]). In our settings, we will use a prime conductor f = q and consider ∆ q = q 2 ∆ K , for a fundamental discriminant ∆ K . In that case, the order of the kernel of φq is given by the classical analytic class number formula (see for instance [START_REF] Buchmann | Binary Quadratic Forms[END_REF])

h(O ∆q ) h(O ∆ K ) = q -(∆ K /q) if ∆ k < -4, (q -(∆ K /q))R ∆ K /R ∆q if ∆ k > 0.
(2)

Note that in the case of real quadratic fields, ǫ ∆q = ǫ t ∆ K for a positive integer t, hence R ∆q /R ∆ K = t and t | (q -(∆ K /q)).

Representation of the Classes

Working with ideals modulo the equivalence relation of the class group is essentially equivalent to work with binary quadratic forms modulo SL 2 (Z) (cf. Section 5.2 of [START_REF] Cohen | A Course in Computational Algebraic Number Theory[END_REF]). Moreover, quadratic forms are more suited to an algorithmic point of view. Every ideal a of O ∆ can be written as a = m aZ

+ -b+ √ ∆ 2 Z with m ∈ Z, a ∈ N and b ∈ Z such that b 2 ≡ ∆ (mod 4a
). In the remainder, we will only consider primitive integral ideals, which are those with m = 1. This notation also represents the binary quadratic form ax 2 + bxy + cy 2 with b 2 -4ac = ∆. This representation of the ideal is unique if the form is normal (see below). We recall here some facts about binary quadratic forms. Definite Positive Forms. Let us first define the crucial notion of reduction. . We denote this form by ρ(f ) and by Rho a corresponding algorithm. The reduction then consists in normalising f , and then iteratively replacing f by ρ(f ) until f is reduced. The time complexity of this (Lagrange-Gauß) algorithm is quadratic (see [START_REF] Buchmann | Binary Quadratic Forms[END_REF]). It returns a reduced form g which is equivalent to f modulo SL 2 (Z). We will call matrix of the reduction, the matrix M such that g = f.M . The reduction procedure yields a uniquely determined reduced form in the class modulo SL 2 (Z). Indefinite Forms. Our main result will deal with forms of positive discriminant. Here is the definition of a reduced indefinite form.

Definition 2. The form f = [a, b, c] is called normal if -a < b ≤ a. It is called reduced if it is normal, a ≤ c,
Definition 3. The form f = [a, b, c] of positive discriminant ∆ is reduced if √ ∆ -2|a| < b < √ ∆ and normal if -|a| < b ≤ |a| for |a| ≥ √ ∆, and √ ∆ -2|a| < b < √ ∆ for |a| < √ ∆.
The reduction process is similar to the definite positive case. The time complexity of the algorithm is still quadratic (see [BV07, Theorem 6.6.4]). It returns a reduced form g which is equivalent to f modulo SL 2 (Z). The main difference with forms of negative discriminant is that there will in general not exist a unique reduced form per class, but several organised in a cycle structure i. e., when f has been reduced then subsequent applications of ρ give other reduced forms. Definition 4. Let f be an indefinite binary quadratic form, the cycle of f is the sequence (ρ i (g)) i∈Z where g is a reduced form which is equivalent to f . From Theorem 6.10.3 from [START_REF] Buchmann | Binary Quadratic Forms[END_REF], the cycle of f consists of all reduced forms in the equivalence class of f . Actually, the complete cycle is obtained by a finite number of application of ρ as the process is periodic. It has been shown in [START_REF] Buchmann | Short Representation of Quadratic Integers[END_REF] that the period length ℓ of the sequence of reduced forms in each class of a class group of discriminant ∆ satisfies R ∆ log ∆ ≤ ℓ ≤ 2R ∆ log 2 + 1. Our factoring algorithm will actually take place in the principal equivalence class. The following definition exhibits the principal form of discriminant ∆.

Definition 5. The reduced form [1, ⌊ √ ∆⌋, (⌊ √ ∆⌋ 2 -∆)/4
] of discriminant ∆ is called the principal form of discriminant ∆, and will be denoted 1 ∆ .

2.3 Reduction of the Forms [q 2 , kq, (k 2 ± p)/4] and Heuristics

In this subsection, p and q are two distinct primes of the same bit-size λ and p ≡ 1 mod 4 (resp. p ≡ 3 mod 4) when we deal with positive (resp. negative) discriminant. Our goal is to factor the numbers pq 2 with the special normalised quadratic forms [q 2 , kq, (k 2 + p)/4] or [q 2 , kq, (k 2 -p)/4], depending whether we work with a negative discriminant ∆ q = -pq 2 or with a positive one ∆ q = pq 2 . If p and q have the same size, these forms are clearly not reduced neither in the imaginary setting nor in real one. But as we shall see, we can find the reduced forms which correspond to the output of the reduction algorithm applied on these forms.

Suppose that we know a form fk , either definite positive or indefinite, which is the reduction of a form f k = [q 2 , kq, (k 2 ± p)/4] where k is an integer. Then fk represents the number q 2 . More precisely, if

M k = α β γ δ ∈ SL 2 (Z) is the matrix such that fk = f k .M k , then fk .M -1 k = f k and q 2 = f k (1, 0) = fk (δ, -γ).
In Section 3, we will see that provided they are relatively small compared to ∆ q , the values δ and -γ can be found in polynomial time with a new variant of Coppersmith method. Our factoring algorithm can be sketched as follows: find such a form fk and if the coefficients of M k are sufficiently small, retrieve δ and -γ and the non-trivial factor q 2 of ∆ q . In this paragraph, we give some heuristics on the size of such a matrix M k and discuss their relevance. If M is a matrix we denote by |M | the max norm, i. e., the maximal coefficient of M in absolute value.

In the imaginary case, it is showed in the proof of [CL09, Theorem 2] that the forms f k belong to different classes of the kernel of the map φq , depending on k, so the reduced equivalent forms fk are the unique reduced elements of the classes of the kernel. To prove the correctness of our attack on NICE, we need the following heuristic1 (indeed, the root finding algorithm of Section 3 recovers roots up to |∆ q | 1/9 ): Heuristic 1 (Imaginary case) Given a reduced element fk of a nontrivial class of ker φq , the matrix of reduction M k is such that |M k | < |∆ q | 1/9 with probability asymptotically close to 1.

From Lemma 5.6.1 of [BV07], |M k | < 2 max{q 2 , (k 2 + p)/4}/ pq 2 . As f k is normalised, |k| ≤ q and |M k | < 2q/ √ p ≈ |∆ q | 1/6
. Note that we cannot reach such a bound with our root finding algorithm. Experimentally, for random k, |M k | can be much smaller. For example, if the bit-size λ of p and q equals 100, the mean value of |M k | is around |∆ q | 1/11.7 . Our heuristic can be explained as follows.

A well-known heuristic in the reduction of positive definite quadratic forms (or equivalently, two-dimensional lattices) is that if [a, b, c] is a reduced quadratic form of discriminant ∆, then a and c should be close to √ ∆. This cannot hold for all reduced forms, but it can be proved to hold for an overwhelming majority of reduced forms. Applied to fk = [a, b, c], this means that we expect a and c to be close to |∆ q | 1/2 . Now, recall that q 2 = fk (δ, -γ) = aδ 2 -bδγ + cγ 2 , which leads to δ and γ close to q 2 /a = q/ √ a ≈ q/|∆ q | 1/4 ≈ |∆ q | 1/12 . Thus, we expect that

|M k | ≤ |∆ q | 1/12
. And this explains why we obtained experimentally the bound |∆ q | 1/11.7 . Figure 1(a) shows a curve obtained by experimentation, which gives the probability that |M k | < |∆ q | 1/9 for random k, in function of λ. This curve also supports our heuristic.

In the real case, we prove in the following theorem that R ∆q /R ∆ K forms f k are principal and we exhibit the generators of the corresponding primitive ideals.

Theorem 1. Let ∆ K be a fundamental positive discriminant, ∆ q = ∆ K q 2 where q is an odd prime conductor. Let ε ∆ K (resp. ε ∆q ) be the fundamental unit of O ∆ K (resp. O ∆q ) and t such that ε t ∆ K = ε ∆q . Then the principal ideals of O ∆q generated by qε i ∆ K correspond to quadratic forms f k(i) = [q 2 , k(i)q, (k(i) 2 -p)/4] with i ∈ {1, . . . , t -1} and k(i) is an integer defined modulo 2q computable from ε i ∆ K mod q.

Proof. Let α i = qε i ∆ K with i ∈ {1, . . . , t -1}. Following the proof of [BTW95, Proposition 2.9], we detail here the computation of a i = α i O ∆q . Let x i and y i be two integers such that ε

i ∆ K = x i + y i ω K . Then α i = qx i + y i q∆ K (1 -q)/2 + y i 1 2 (∆ q + ∆ q ),
and α i is an element of O ∆q . Let m i , a i and b i be three integers such that

a i = m i a i Z + -bi+ √ ∆q 2
. As mentioned in the proof of [BTW95, Proposition 2.9], m i is the smallest positive coefficient of ∆ q /2 in a i . As O ∆q is equal to Z + (∆ q + ∆ q )/2Z, α i O ∆q is generated by α i and α i (∆ q + ∆ q )/2 as a Z-module. So a simple calculation gives that m i = gcd(y i , q(x i + y i ∆ K /2)). As ε i ∆ K is not an element of O ∆q , we have gcd(y i , q) = 1 so m i = gcd(y i , x i + y i ∆ K /2). The same calculation to find m ′ i for the ideal

ε i ∆ K O ∆ K reveals that m i = m ′ i . As ε i ∆ K O ∆ K = O ∆ K we must have m ′ i = 1. Now, N(a i ) = |N(α i )| = q 2 and N(a i ) = m 2
i a i = a i and therefore a i = q 2 . Let us now find b i . Note that b i is defined modulo 2a i . Since α i ∈ α i O ∆q , there exist µ i and ν i such that α i = a i µ i + (-b i + ∆ q )/2ν i . By identification in the basis (1, ∆ q ), ν k = 1 and by a multiplication by 2, we obtain 2qx i + qy i ∆ K ≡ -b i y i (mod 2a i ). As b i ≡ ∆ q (mod 2), we only have to determine b i modulo q 2 . As y i is prime to q, we have b i ≡ k(i)q (mod q 2 ) with k(i) ≡ -2x i /y i -∆ K (mod q). Finally, as we must have -a i < b ≤ a i if a i > ∆ q and else ∆ q -2a i < b < ∆ q , k(i) is the unique integer with k(i) ≡ ∆ q (mod 2) and k(i) ≡ -2x i /y i -∆ K (mod q), such that b = k(i)q satisfies that inequalities. Eventually, the principal ideal of O ∆q generated by qε i ∆ K corresponds to the form [q 2 , k(i)q, c i ] with

c i = (b 2 i -∆ q )/(4a i ) = (k(i) 2 -∆ K )/4. ⊓ ⊔
From this theorem, we see that if we go across the cycle of principal forms, then we will find reduced forms fk . To analyse the complexity of our factoring algorithm, we have to know the distribution of these forms on the cycle. An appropriate tool is the Shanks distance d (see [BV07, Definition 10.1.4]) which is close to the number of iterations of Rho between two forms. One has d(1 ∆q , f k(i) ) = iR ∆ K . From Lemma 10.1.8 of [START_REF] Buchmann | Binary Quadratic Forms[END_REF], |d( fk(i) , f k(i) )| < log q, for all i = 1, 2, . . . , t -1. Let j be the smallest integer such that 0 < jR ∆ K -2 log q, then as jR

∆ K = d(f k(i) , f k(i+j) ) = d(f k(i) , fk(i) ) + d( fk(i) , fk(i+j) ) + d( fk(i+j) , f k(i+j)
), from the triangle inequality, one has jR ∆ K < 2 log(q) + |d( fk(i) , fk(i+j) )|. So, |d( fk(i) , fk(i+j) )| > jR ∆ K -2 log q > 0. This inequality proves that f k(i) and f k(i+j) do not reduce to the same form. Experiments actually show that asymptotically, |d( fk(i) , f k(i) )| is very small on average (smaller than 1). As a consequence, as pictured in figure 2,d

(1 ∆q , fk(i) ) ≈ iR ∆ K . R∆ K f k(1) fk(1) f k(2) fk(2) f k(3) fk(3) 1∆ q

Fig. 2. Repartition of the forms fk(i) along the principal cycle

Moreover, as in the imaginary case, experiments show that asymptotically the probability that the norm of the matrices of reduction, |M k | is smaller than ∆ 1/9 q is close to 1 (see figure 1(b)). This leads to the following heuristic.

Heuristic 2 (Real case) From the principal form 1 ∆q , a reduced form fk such that the matrix of the reduction, M k , satisfy

|M k | < ∆ 1/9 q , can be found in O(R ∆ K ) successive applications of Rho.
We did also some experiments to investigate the case where the bit-sizes of p and q are unbalanced. In particular when the size of q grows, the norm of the matrix of reduction becomes larger. For example, for a 100-bit p and a 200-bit q (resp. a 300-bit q), more than 95% (resp. 90%) of the fk have a matrix

M k with |M k | < ∆ 1/6.25 q (resp. |M k | < ∆ 1/5.44 q ).

A Rigorous Homogeneous Variant of Coppersmith's Root Finding Method

Our factoring algorithm searches many times for small modular roots of degree two homogeneous polynomials and the most popular technique to find them is based on Coppersmith's method (see [START_REF] Coppersmith | Small Solutions to Polynomial Equations, and Low Exponent RSA Vulnerabilities[END_REF] or May's survey [START_REF] May | Using LLL-Reduction for Solving RSA and Factorization Problems: A Survey[END_REF]).

Our problem is the following: Given f (x, y) = x 2 + bxy + cy 2 a (monic) binary quadratic form and N = pq 2 an integer of unknown factorisation, find (x 0 , y 0 ) ∈ Z 2 such that f (x 0 , y 0 ) ≡ 0 (mod q 2 ), while |x 0 |, |y 0 | ≤ M , where M ∈ N. The usual technique for this kind of problems is only heuristic, since it is the gcd extension of bivariate congruences. Moreover, precise bounds cannot be found in the litterature. Fortunately, because our polynomial is homogeneous, we will actually be able to prove the method. This homogenous variant is quite similar to the one-variable standard Coppersmith method, but is indeed even simpler to describe and more efficient since there is no need to balance coefficients. We denote as • the usual Euclidean norm for polynomials. The main tool to solve this problem is given by the following variant of the widespread elementary Howgrave-Graham's lemma [START_REF] Howgrave-Graham | Finding small roots of univariate modular equations revisited[END_REF].

Lemma 1. Let B ∈ N and g(x, y) ∈ Z[x, y] be a homogeneous polynomial of total degree δ. Let M > 0 be a real number and suppose that ||g(x, y)|| <

B √ δ+1M δ
then for all x 0 , y 0 ∈ Z such that g(x 0 , y 0 ) ≡ 0 (mod B) and |x 0 |, |y 0 | ≤ M , g(x 0 , y 0 ) = 0.

Proof. Let g(x, y) = δ i=0 g i x i y δ-i where some g i s might be zero. We have

|g(x 0 , y 0 )| ≤ δ i=0 |g i ||x i 0 y δ-i 0 | ≤ M δ δ i=0 |g i | ≤ M δ √ δ + 1 g(x, y) < B
and therefore g(x 0 , y 0 ) = 0.

⊓ ⊔

The trick is then to find only one small enough bivariate homogeneous polynomial satisfying the conditions of this lemma and to extract the rational root of the corresponding univariate polynomial with standard techniques. On the contrary, the original Howgrave-Graham's lemma suggests to look for two polynomials of small norm having (x 0 , y 0 ) as integral root, and to recover it via elimination theory. The usual way to obtain these polynomials is to form a lattice spanned by a special family of polynomials, and to use the LLL algorithm (cf. [START_REF] Lenstra | Factoring Polynomials with Rational Coefficients[END_REF]) to obtain the two "small" polynomials. Unfortunately, this For the case we are most interested in, δ = 2, N = pq 2 with p and q of the same size, i. e., α = 2/3 then λ = 3/2 and we can asymptotically get roots up to N β with β = 1 9 . If we take m = 4 and t = 6, i. e., we work with a lattice of dimension 13, we get from (3) that β ≈ 1 10.63 and with a 31-dimensional lattice (m = 10 and t = 15), β ≈ 1 9.62 . If the size of q grows compared to p, i. e., α increases towards 1, then β increases towards 1/4. For example, if q is two times larger than p, i. e., α = 4/5 then β = 1/6.25. For α = 6/7, we get β ≈ 1/5.44.

We will call HomogeneousCoppersmith the algorithm which implements this method. It takes as input an integer N = pq 2 and a binary quadratic form [a, b, c], from which we deduce the unitary polynomial x 2 +b ′ xy +c ′ y 2 , by dividing both b and c by a modulo N , and the parameters m and t. In fact, this method will only disclose proper representations of q 2 , those for which x and y are coprime, but we note that f k properly represents q 2 , and therefore so does our form [a, b, c].

The case α = 1 of Theorem 2 can already be found in Joux's book [START_REF] Joux | Algorithmic Cryptanalysis[END_REF] and we mention that a similar technique has already been independently investigated by Bernstein in [START_REF] Bernstein | List decoding for binary Goppa codes[END_REF].

4 A Õ(p 1/2 )-Deterministic Factoring Algorithm for pq 2
We detail our new quadratic form-based factoring algorithm for numbers of the form pq 2 . In this section, p and q will be of same bit-size, and p ≡ 1 (mod 4).

The Algorithm

Roughly speaking, if ∆ q = N = pq 2 , our factoring algorithm, depicted in Fig. 3, exploits the fact that the non-reduced forms f k = [q 2 , kq, -] reduce to forms fk for which there exists a small pair (x 0 , y 0 ) such that q 2 | fk (x 0 , y 0 ) while q 2 | N . From Theorem 1, we know that these reduced forms appear on the principal cycle of the class group of discriminant ∆ q . To detect them, we start a walk in the principal cycle from the principal form 1 N , and apply Rho until the Coppersmith-like method finds these small solutions.

Input: N = pq 2 , m, t Output: p, q 1. h ← 1N 2. while (x0, y0) not found do 2.1. h ← Rho(h) 2.
2. x0/y0 ← HomogeneousCoppersmith(h, N, m, t) 3. q ← Sqrt(Gcd(h(x0, y0), N )) 4. return (N/q 2 , q) Fig. 3. Factoring N = pq 2

Heuristic Correctness and Analysis of Our Algorithm

Assuming Heuristic 2, starting from 1 N , after O(R p ) iterations, the algorithm will stop on a reduced form whose roots will be found with our Coppersmithlike method (for suitable values of m and t) since they will satisfy the expected N 1/9 bound. The computation of gcd(h(x 0 , y 0 ), N ) will therefore expose q 2 and factor N . The time complexity of our algorithm is then heuristically O(R p Poly(log N )), whereas the space complexity is O(log N ). The worst-case complexity is O(p 1/2 log p Poly(log N )). For small regulators, such as in REAL-NICE cryptosystem (see. Subsection 5.1), the time complexity is polynomial.

This algorithm can be generalised with a few modifications to primes p such that p ≡ 3 (mod 4), by considering ∆ q = 4pq 2 . Moreover if the bit-sizes of p and q are unbalanced, our experiments suggest that the size of the roots will be small enough (see end of Subsection 2.3 and Section 3), so the factoring algorithm will also work in this case, with the same complexity.

Comparison with other Deterministic Factorisation Methods. Boneh, Durfee and Howgrave-Graham presented in [START_REF] Boneh | Factoring N = p r q for large r[END_REF] an algorithm for factoring integers N = p r q. Their main result is the following: Lemma 3 ([BDH99]). Let N = p r q be given, and assume q < p c for some c. Furthermore, assume that P is an integer satisfying |P -p| < p 1-c r+c -2 r d . Then the factor p may be computed from N , r, c and P by an algorithm whose running time is dominated by the time it takes to run LLL on a lattice of dimension d.

For r = 2 and c = 1, this leads to a deterministic factoring algorithm which consists in exhaustively search for an approximation P of p and to solve the polynomial equation (P + X) 2 ≡ 0 (mod p 2 ) with a method à la Coppersmith. The approximation will be found after O(p 1/3 ) = O(N 1/9 ) iterations.

The fastest deterministic generic integer factorisation algorithm is actually a version of Strassen's algorithm [START_REF] Strassen | Einige Resultate über Berechnungskomplexität[END_REF] from Bostan, Gaudry and Schost [START_REF] Bostan | Linear Recurrences with Polynomial Coefficients and Application to Integer Factorization and Cartier-Manin Operator[END_REF], who ameliorates a work of Chudnovsky and Chudnovsky [START_REF] Chudnovsky | Approximations and Complex Multiplication According to Ramanujan[END_REF] and proves a complexity of O(M int ( 4 √ N log N )) where M int is a function such that integers of bit-size d can me multiplied in M int (d) bit operations. More precisely, for numbers of our interest, Lemma 13 from [START_REF] Bostan | Linear Recurrences with Polynomial Coefficients and Application to Integer Factorization and Cartier-Manin Operator[END_REF] gives the precise complexity: In particular, for b = N 1/3 , the complexity is Õ(N 1/6 ), with a very large space complexity compared to our algorithm. Moreover, none of these two last of algorithms can actually factor an integer of cryptographic size. The fact that a prime divisor has a small regulator does not help in these algorithms, whereas it makes the factorisation polynomial in our method.

Cryptanalysis of the NICE Cryptosystems

Hartmann, Paulus and Takagi proposed the elegant NICE encryption scheme (see [HPT99,PT99,PT00]), based on imaginary quadratic fields and whose main feature was a quadratic decryption time. Later on, several other schemes, including (special) signature schemes relying on this framework have been proposed. The public key of these NICE cryptosystems contains a discriminant ∆ q = -pq 2 together with a reduced ideal h whose class belongs to the kernel of φq . The idea underlying the NICE cryptosystem is to hide the message behind a random element [h] r of the kernel. Applying φq will make this random element disappear, and the message will then be recovered.

In [START_REF] Jacobson | An Adaptation of the NICE Cryptosystem to Real Quadratic Orders[END_REF], Jacobson, Scheidler and Weimer embedded the original NICE cryptosystem in real quadratic fields. Whereas the idea remains essentially the same as the original, the implementation is very different. The discriminant is now ∆ q = pq 2 , but because of the differences between imaginary and real setting, these discriminant will have to be chosen carefully. Among these differences, the class numbers are expected to be small with very high probability (see the Cohen-Lenstra heuristics [START_REF] Cohen | Heuristics on class groups[END_REF]). Moreover, an equivalence class does not contain a unique reduced element anymore, but a multitude of them, whose number is governed by the size of the fundamental unit. The rough ideas to understand these systems and our new attacks are given in the following. The full description of the systems is omitted for lack of space but can be found in [HPT99, JSW08].

Polynomial-Time Key Recovery in the Real Setting

The core of the design of the REAL-NICE encryption scheme is the very particular choice of the secret prime numbers p and q such that ∆ K = p and ∆ q = pq 2 . They are chosen such that the ratio R ∆q /R ∆ K is of order of magnitude of q and that R ∆ K is bounded by a polynomial in log(∆ K ). To ensure the first property, it is sufficient to choose q such that q -∆ K q is a small multiple of a large prime. If the second property is very unlikely to naturally happen since the regulator of p is generally of the order of magnitude of √ p, it is indeed quite easy to construct fundamental primes with small regulator. The authors of [START_REF] Jacobson | An Adaptation of the NICE Cryptosystem to Real Quadratic Orders[END_REF] suggest to produce a prime p as a so-called Schinzel sleeper, which is a positive squarefree integer of the form p = a 2 x 2 + 2bx + c with a, b, c, x in Z, a = 0 and b 2 -4ac dividing 4 gcd(a 2 , b) 2 . Schinzel sleepers are known to have a regulator of the order log(p) (see [START_REF] Cheng | Some Results Concerning Certain Periodic Continued Fractions[END_REF]). Some care must be taken when setting the (secret) a, b, c, x values, otherwise the resulting ∆ q = pq 2 is subject to factorisation attacks described in [START_REF] Weimer | An Adaptation of the NICE Cryptosystem to Real Quadratic Orders[END_REF]. We do not provide here more details on these choices since the crucial property for our attack is the fact that the regulator is actually of the order log(p). The public key consists of the sole discriminant ∆ q . The message is carefully embedded (and padded) into a primitive O ∆q -ideal so that it will be recognised during decryption. Instead of moving the message ideal m to a different equivalence class (like in the imaginary case), the encryption actually hides the message in the cycle of reduced ideal of its own equivalent class by multiplication of a random principal O ∆q -ideal h (computed during encryption). The decryption process consists then in applying the (secret) map φq and perform an exhaustive search for the padded message in the small cycle of φq ([mh]). This exhaustive search is actually possible thanks to the choice of p which has a very small regulator. Like in the imaginary case, the decryption procedure has a quadratic complexity and significantly outperforms an RSA decryption for any given security level (see Table 3 from [START_REF] Jacobson | An Adaptation of the NICE Cryptosystem to Real Quadratic Orders[END_REF]).

Unfortunately, due to the particular but necessary choice of the secret prime p, the following result states the total insecurity of the REAL-NICE system.

Result 1 Algorithm 3 recovers the secret key of REAL-NICE in polynomial time in the security parameter under Heuristic 2 since the secret fundamental discriminant p is chosen to have a regulator bounded by a polynomial in log p.

We apply the cryptanalysis on the following example. The Schinzel polynomial S(X) = 2725 2 X 2 + 2 • 3815X + 2 produces a suitable 256-bit prime p for the value X 0 = 103042745825387139695432123167592199. This prime has a regulator R ∆ K ≃ 90.83. The second 256-bit prime q is chosen following the recommendations from [START_REF] Weimer | An Adaptation of the NICE Cryptosystem to Real Quadratic Orders[END_REF]. This leads to a the discriminant ∆q = 28736938823310044873380716142282073396186843906757463274792638734144060602830510 80738669163489273592599054529442271053869832485363682341892124500678400322719842 63278692833860326257638544601057379571931906787755152745236263303465093

Our algorithm recovers the prime q = 60372105471499634417192859173853663456123015267207769653235558092781188395563 from ∆ q after 45 iterations in 42.42 seconds on a standard laptop. The rational root is x0 y0 equal to -2155511611710996445623 3544874277134778658948 , where x 0 and y 0 satisfy log(∆q) log(|x0|) ≃ 10.8 and log(∆q) log(|y0|) ≃ 10.7.

Polynomial-Time Key Recovery of the Original NICE

As mentioned above, the public key of the original NICE cryptosystem contains the representation of a reduced ideal h whose class belongs to the kernel of the surjection φq . The total-break of the NICE cryptosystem is equivalent to solving the following kernel problem.

Definition 6 (Kernel Problem [START_REF] Biehl | Efficient Undeniable Signature Schemes based on Ideal Arithmetic in Quadratic Orders[END_REF]). Let λ be an integer, p and q be two λ-bit primes with p ≡ 3 (mod 4). Fix a non-fundamental discriminant ∆ q = -pq 2 . Given an element [h] of ker φq , factor the discriminant ∆ q .

Castagnos and Laguillaumie proposed in [START_REF] Castagnos | On the Security of Cryptosystems with Quadratic Decryption: The Nicest Cryptanalysis[END_REF] a polynomial-time algorithm to solve this problem. We propose here a completely different solution within the spirit of our factorisation method and whose complexity is also polynomial-time. As discuss in Subsection 2.3, the idea is to benefit from the fact that the public ideal h corresponds to a reduced quadratic form, fk , which represents q 2 . We thus find these x 0 and y 0 such that gcd( fk (x 0 , y 0 ), ∆ q ) = q 2 with the Coppersmith method of Section 3.

Definition 1 .

 1 A binary quadratic form f is a degree 2 homogeneous polynomial f (x, y) = ax 2 +bxy +cy 2 where a, b and c are integers, and is denoted by [a, b, c]. The discriminant of the form is ∆ = b 2 -4ac. If a > 0 and ∆ < 0, the form is called definite positive and indefinite if ∆ > 0. Let M ∈ SL 2 (Z) with M = α β γ δ , and f = [a, b, c], a binary quadratic form, then f.M is the equivalent binary quadratic form f (αx + βy, γx + δy).

  and if b ≥ 0 for a = c. The procedure which transforms a form f = [a, b, c] into a normal one consists in setting s such that b + 2sa belongs to the right interval (see [BV07, (5.4)]) and producing the form [a, b+2sa, as 2 +bs+c]. Once a form f = [a, b, c] is normalised, a reduction step consists in normalising the form [c, -b, a]

Fig. 1 .

 1 Fig. 1. Probability that |M k | < |∆q| 1/9 in function of the bit-size λ of p and q

  Lemma 4 ([BGS07]). Let b, N be two integers with 2 ≤ b < N . One can compute a prime divisor of N bounded by b, or prove that no such divisor exists in O M int ( √ b log N ) + log bM int (log N ) log log N bit operations and space O( √ b log N ) bits.

In the full version, we prove a probabilistic version of Heuristic 1.

Acknowledgements. We warmly thank Denis Simon and Brigitte Vallée for helpful discussions and the reviewers for their useful comments. Part of this work was supported by the Commission of the European Communities through the ICT program under contract ICT-2007-216676 ECRYPT II.

⋆ This work was done while

reduction does not guarantee that these polynomials will be algebraically independent, and the elimination can then lead to a trivial relation. Consequently, this bivariate approach is heuristic. Fortunately, for homogeneous polynomials, we can take another approach by using Lemma 1 and then considering a univariate polynomial with a rational root. This makes the method rigorous and slightly simpler since we need a bound on g(x, y) and not on g(xX, yY ) if X and Y are bounds on the roots and therefore the resulting lattice has smaller determinant than in the classical bivariate approach.

To evaluate the maximum of the bound we can obtain, we need the size of the first vector provided by LLL which is given by: Lemma 2 (LLL). Let L be a full-rank lattice in Z d spanned by an integer basis B = {b 1 , . . . , b d }. The LLL algorithm, given B as input, will output a non-zero vector u ∈ L satisfying u ≤ 2 (d-1)/4 det(L) 1/d in time O(d 6 log 3 (max b i )).

We will now prove the following general result regarding the modular roots of bivariate homogeneous polynomials which can be of independent interest. Theorem 2. Let f (x, y) ∈ Z[x, y] be a homogeneous polynomial of degree δ with f (x, 0) = x δ , N be a non-zero integer and α be a rational number in [0, 1], then one can retrieve in polynomial time in log N , δ and the bit-size of α, all the rationals x 0 /y 0 , where x 0 and y 0 are integers such that gcd(f (x 0 , y 0 ),

Proof. Let b be a divisor of N for which their exists (x 0 , y 0 ) ∈ Z 2 such that b = gcd(f (x 0 , y 0 ), N ) ≥ N α . We define some integral parameters (to be specified later) m, t and t ′ with t = m + t ′ and construct a family of δt + 1 homogeneous polynomials g and h of degree δt such that (x 0 , y 0 ) is a common root modulo b m . More precisely, we consider the following polynomials g i,j (x, y) = x j y δ(t-i)-j f i N m-i for i = 0, . . . , m -1, j = 0, . . . , δ -1 h i (x, y) = x i y δt ′ -i f m for i = 0, . . . , δt ′ .

We build the triangular matrix L of dimension δt + 1, containing the coefficients of the polynomials g i,j and h i . We will apply LLL to the lattice spanned by the rows of L. The columns correspond to the coefficients of the monomials y δt , xy δt-1 , . . . , x δt-1 y, x δt . Let β ∈ [0, 1] such that M = N β . The product of the diagonal elements gives det(L) = N δm(m+1)/2 . If we omit the quantities that do not depend on N , to satisfy the inequality of Lemma 1 with the root bound M , the LLL bound from Lemma 2 implies that we must have δm(m + 1)/2 ≤ (δt + 1)(αm -δtβ)

(3) and if we set λ such that t = λm, this gives asymptotically β ≤ α δλ -1 2δλ 2 , which is maximal when λ = 1 α , and in this case, β max = α 2 /(2δ). The vector output by LLL gives a homogeneous polynomial f (x, y) such that f (x 0 , y 0 ) = 0 thanks to Lemma 1. Let r = x/y, any rational root of the form x 0 /y 0 can be found by extracting the rational roots of f ′ (r) = 1/y δt f (x, y) with classical methods. ⊓ ⊔

Result 2 The Homogeneous Coppersmith method from Section 3 solves the Kernel Problem in polynomial time in the security parameter under Heuristic 1.

We apply our key recovery on the example of NICE proposed in [JJ00,CL09]: The public key consists in ∆ q and h = (a, b). Our Coppersmith method finds in less that half a second the root u 0 = -103023911 349555951 = x0 y0 and h(x0, y0) = 5363123171977038839829609999282338450991746328236957351089 4245774887056120365979002534633233830227721465513935614971 593907712680952249981870640736401120729 = q 2 .

All our experiments have been run on a standard laptop under Linux with software Sage. The lattice reduction have been performed with Stehlé's fplll [Ste].


Yongchao Xu

Edwin Carlinet

Thierry Géraud

Laurent Najman

MEANINGFUL DISJOINT LEVEL LINES SELECTION

Keywords: Mathematical morphology, Level lines, Tree of shapes, Saliency map, Connected filters

c: Bottom-up level lines selection; 704 level lines. d: Top-down selection; 404 level lines. e: Selection by gradient's magnitude; 181 level lines. Fig. 1: An example of the proposed method. (b-d): selected level lines.

INTRODUCTION

Many authors [START_REF] Caselles | Topographic maps and local contrast changes in natural images[END_REF][START_REF] Cao | Extracting meaningful curves from images[END_REF][START_REF] Lisani | On the theory of planar shape[END_REF] claim that significant contours of objects coincide with segments of level lines in images. Each level line is the contour of a level set or a shape which is a connected component with holes filled. The set of shapes equipped with the inclusion relationship allows to represent an image by a tree structure, called tree of shapes [START_REF] Monasse | Fast computation of a contrast-invariant image representation[END_REF] or topographic map [START_REF] Caselles | Topographic maps and local contrast changes in natural images[END_REF]. Numerous applications using the tree of shapes have been proposed: image filtering (grain filter [START_REF] Caselles | Grain filters[END_REF], our proposed morphological shaping [START_REF] Xu | Morphological Filtering in Shape Spaces: Applications using Tree-Based Image Representations[END_REF]), image simplification [START_REF] Cao | Extracting meaningful curves from images[END_REF][START_REF] Ballester | Level lines selection with variational models for segmentation and encoding[END_REF][START_REF] Xu | Salient level lines selection using the mumford-shah functional[END_REF], image segmentation [START_REF] Pardo | Semantic image segmentation using morphological tools[END_REF][START_REF] Cardelino | Region based segmentation using the tree of shapes[END_REF][START_REF] Xu | Context-based energy estimator: Application to object segmentation on the tree of shapes[END_REF], scenery images analysis [START_REF] Song | Locating image background by monotonic tree[END_REF][START_REF] Song | Analyzing scenery images by monotonic tree[END_REF], object recognition [START_REF] Pan | Preferential image segmentation using trees of shapes[END_REF], and texture indexing [START_REF] Xia | Shape-based invariant texture indexing[END_REF]. In spite of all these abundant applications, the strong potential of the tree of shapes is still under-exploited.

The number of connected components in the tree of shapes is about as large as the number of pixels in the image. Hence, it is not straightforward to know a priori the structure of the tree (as shown in Fig. 1a). However, this information is fundamental for a deep tree analysis and developing tree-based methods. In this paper, we propose an efficient algorithm that selects a set of disjoint level lines of the tree of shapes from meaningful ones to meaningless ones. These level lines represent the main structure of the original tree, which can now be easily perceived from the contour clean image reconstructed from those level lines. This is the main contribution of this paper. An example of selected level lines is given in Fig. 1: In (b), too many level lines are preserved that makes it difficult to perceive the main structure of the image. In (c) and (d), the bottom-up and top-down level lines selection do not retrieve correctly the image geometry, whereas the selection by meaningfulness in (e) reveals the main structure of the image.

There exist many similar applications relying on the ex- traction of level lines from the tree of shapes to achieve image filtering, simplification, and segmentation tasks. Salembier and Wilkinson [START_REF] Salembier | Connected operators[END_REF] give a detailed review of the tree-based filtering strategies. The works in [START_REF] Pardo | Semantic image segmentation using morphological tools[END_REF][START_REF] Cardelino | Region based segmentation using the tree of shapes[END_REF][START_REF] Desolneux | Edge detection by helmholtz principle[END_REF][START_REF] Cao | Extracting meaningful curves from images[END_REF][START_REF] Xu | Context-based energy estimator: Application to object segmentation on the tree of shapes[END_REF][START_REF] Ballester | Level lines selection with variational models for segmentation and encoding[END_REF][START_REF] Pan | Top-down image segmentation using the Mumford-Shah functional and level set image representation[END_REF][START_REF] Xu | Salient level lines selection using the mumford-shah functional[END_REF] are closest to what we propose here. In [START_REF] Pardo | Semantic image segmentation using morphological tools[END_REF], the author proposes a segmentation algorithm that selects the perceptible level lines matching some criteria, e.g., number of T-junctions, compactness, and contrast. In [START_REF] Cardelino | Region based segmentation using the tree of shapes[END_REF], the authors remove the level lines that enclose a region similar to its parent w.r.t. a histogrambased distance. Then they select interesting regions by identifying parts of the tree having an homogeneous histogram. Desolneux et al. [START_REF] Desolneux | Edge detection by helmholtz principle[END_REF] and Cao et al. [START_REF] Cao | Extracting meaningful curves from images[END_REF] define the meaningfulness of a given level line using the a contrario model based on the computation of a number of false alarms (NFA). In [START_REF] Xu | Context-based energy estimator: Application to object segmentation on the tree of shapes[END_REF], we also propose a context-based energy estimator to measure the significance. In [START_REF] Ballester | Level lines selection with variational models for segmentation and encoding[END_REF][START_REF] Pan | Top-down image segmentation using the Mumford-Shah functional and level set image representation[END_REF][START_REF] Xu | Salient level lines selection using the mumford-shah functional[END_REF], the authors proposed to extract level lines by minimizing piecewise-constant Mumford-Shah functional. To the best of our knowledge, there is no method in the literature that tries to extract a subset of level lines for visualizing the tree. The most similar one might be the tree of bilinear level lines in the chapter 7 of [START_REF] Caselles | Geometric Description of Images as Topographic Maps[END_REF]. The authors use a bilinear interpolation to compute a continuous version of the image. So level lines from different gray levels do not touch each other, which permits a multiscale visualization. The rest of the paper is organized as follows. Some background information about the tree of shapes and Géraud et al. [START_REF] Géraud | A quasi-linear algorithm to compute the tree of shapes of nD images[END_REF]'s tree computation algorithm are provided in Section 2. Our method is detailed in Section 3, and we present in Section 4 some experimental results. Finally, we conclude and give some perspectives in Section 5.

D E B A C F O (a) Input image. A O F B C D E > 2 > 2 > 2 > 2 > 0 > 2 > 1 > 4 (b) Tree of shapes. B A, F C D E O 3 2 0 1 (c) Max-tree.

BACKGROUND

The Tree of Shapes

Let f be an image defined on a domain Ω and with values on an ordered set V (typically R or Z). For any λ ∈ V , the upper level sets X λ and lower level sets X λ of an image f are respectively defined by lower level sets have a natural inclusion structure: ∀ λ 1 ≤ λ 2 , X λ1 ⊇ X λ2 and X λ1 ⊆ X λ2 , which leads to two distinct and dual representations of an image, the Max-tree and the Min-tree [START_REF] Salembier | Antiextensive connected operators for image and sequence processing[END_REF]. The tree of shapes is a fusion of the Maxtree and Min-tree via the notion of shapes [START_REF] Monasse | Fast computation of a contrast-invariant image representation[END_REF]. A shape is defined as a connected component of an upper or lower level set with its holes filled in. Thanks to the inclusion relationship of both kinds of level sets, the set of shapes can be structured into a tree structure, called the tree of shapes. This tree features several interesting properties: it is invariant to contrast changes and forms a self-dual, non-redundant, and complete representation of an image. Furthermore, such a tree inherently embeds a morphological scale-space (the parent of a node/shape is a larger shape) w.r.t. the theory of multi-scale analysis [START_REF] Guigues | Scalesets image analysis[END_REF]. An example of these trees is depicted in Fig. 2.

X λ (f ) = {p ∈ Ω | f (p) ≥ λ} and X λ (f ) = {p ∈ Ω | f (p) ≤ λ}. Both upper and

Tree computation and representation

Though many tree of shapes computation algorithms exist [START_REF] Monasse | Fast computation of a contrast-invariant image representation[END_REF][START_REF] Song | A topdown algorithm for computation of level line trees[END_REF][START_REF] Caselles | Geometric Description of Images as Topographic Maps[END_REF], we rely on the one proposed in [START_REF] Géraud | A quasi-linear algorithm to compute the tree of shapes of nD images[END_REF] that ensures a worst-case quasi-linear complexity. Even if this paper does not aim at explaining this algorithm since it has already been fully described by its authors, we still need to remind that it does not work directly on the domain Ω but on a simplicial version of the 2D discrete grid: the Khalimsky grid. We note K Ω , the domain Ω immersed on this grid. In Fig. 3, the original points of the image are the 2-faces, the boundaries are materialized with 0-faces and 1-faces. Géraud at al. [START_REF] Géraud | A quasi-linear algorithm to compute the tree of shapes of nD images[END_REF]'s algorithm ensures that shapes are open connected sets (e.g., the purple shape in Fig. 3) and that levels lines (shape's borders) are composed of 0-faces and 1-faces only (e.g., the dark curve in Fig. 3). Given a connected open set A, we note Ā the closure of A and ∂A the frontier of Ā, i.e., the level line associated to the shape A.

Using the tree representation from [START_REF] Géraud | A quasi-linear algorithm to compute the tree of shapes of nD images[END_REF], a node is represented by a single pixel (2-face) called the canonical element. The tree is encoded through an image parent : K Ω → K Ω that states the parenthood relationship between nodes. In parent, each non-canonical element point is attached to the canonical element representing the node it belongs to. In the following, we denote by getCanonical : K Ω → K Ω , the routine that returns the canonical element of each point in the image. We also rely on a sorted vector of pixels S such that processing S in the direct order is a top-down traversal of the tree and the reverse order is a bottom-up traversal.

DISJOINT LEVEL LINES SELECTION

In natural images, there are about as many level lines as pixels, which make the tree of shapes difficult to be visualized for further analysis. In fact, most shapes only differ from a few pixels when climbing the tree of shapes, thus many level lines overlap. In this section, we propose an efficient algorithm that allows to select a set of disjoint level lines whatever the selection ordering is.

Disjoint nodes computation

As shortly reviewed in Section 2.2, a shape A is an open connected set in K Ω and the corresponding level line is the border ∂A composed of 0 and 1-faces only. In this section, we focus on an algorithm that enables to access for each shape, the shapes having part of borders in common. In other words, for any shape A, we aim at computing the set of shapes B such that for each B ∈ B, either A ⊂ B ∧ ∂A ∩ ∂B = ∅ or B ⊂ A ∧ ∂A ∩ ∂B = ∅. Actually, these relations are dual. We can focus on computing only the first one: the set of shapes in the parenthood sharing a border with ∂A. Let SES(A) denotes the Smallest Enclosing Shape of Ā. Then, SES(A) is the smallest open set that totally includes the border of A. Since SES(A) is open, its border does not belong to SES(A) and thus, ∂SES(A) ∩ ∂A = ∅. In other words, the level line ∂SES(A) is the first (w.r.t the inclusion relationship) disjoint level line of ∂A. On the contrary, for each B ∈ [A SES(A)[ (where x y is the path from x to y in the tree), the borders of A and B intersect.

Computing SES : K Ω → K Ω for each shape is a standard incremental attribute computation. This is given in Algorithm 1. For each 2-face x in K Ω , we retain the lowest node which contains x (i.e., the closure of x). This is straightforward if the depth of each point in the tree is given. When processing each 0 and 1-face e connecting to x, we have two cases: either e is an internal edge or an external edge of the shape containing x. In the former case, we have depth(e) ≥ depth(x), thus e is ignored. In the latter case, we have depth(e) < depth(x). The node that contains e is a potential smallest enclosing shape. Since Ā = ∪ {x | x ∈ A}, we propagate the SES from the bottom to the top of the tree incrementally, and eventually we get the SES for each node.

Final disjoint nodes selection

Based on the Smallest Enclosing Shapes, we are able to select a set of nodes in an arbitrary order such that any pair of nodes have disjoint borders. In this paper, we consider three selection strategies: top-down, bottom-up and another one based on a meaningfulness criterion. Let Order denotes the sequence of shapes sorted by one of the previous strategy. Order is actually an array of canonical elements. The final disjoint nodes selection is given in Algorithm 2 which proceeds as follows. Every shape has a three-state status: Null, the shape has not yet been processed, it is a candidate shape; Active, the shape has been selected; Unactive, the shape has been dismissed. At first, every shape has the Null status and we process the candidate shapes in the order given by Order. Given a candidate shape A, we have to check that the border does not intersect with an already Active shape, that is: A has not been marked Unactive by a sub-shape and that any shape in [A SES(A)[ is not Active. If both conditions are fulfilled, A is selected, marked Active and we dismiss every shape in ]A SES(A)[ since they intersect A's border. The smallest enclosing shapes computation described in Algorithm 1 has a linear complexity. For the disjoint level lines selection shown in Algorithm 2, the bottom-up selection has also a linear complexity, whereas the other two strategies have a worst case O(N 2 ) complexity, where N is the number of shapes. We have implemented the proposed method using our C++ image processing library [START_REF] Levillain | Why and how to design a generic and efficient image processing framework: The case of the Milena library[END_REF]. Processing a 640 × 960 pixels image takes about 3s on a regular PC station, where the tree computation takes half of the time.

EXPERIMENTAL RESULTS

We have experienced three different orders for disjoint level lines selection: top-down (from the root to the leaves), bottom-up (from the leaves to the root), and by meaningfulness, e.g., the average of the gradient's magnitude along the level lines. An example of such disjoint level lines selec- tions is illustrated in Fig. 4 and Fig. 1. Usually, the selection order based on meaningfulness is preferred, since it preserves better the meaningful level lines. Consequently, it offers better understanding of the tree structure for further analysis.

In Fig. 5, we show an illustration of this method on a bee board image that displays disjoint black hexagons representing places without honey. The simplified image reconstructed from the extracted disjoint level lines is shown in Fig. 5b. The selection is based on the decreasing order of average of the gradient's magnitude along the level lines. The main structure of the tree of shapes can be easily perceived in Fig. 5b and Fig. 5d. By removing those small nodes, we obtain the object segmentation result presented in Fig. 5e. In Fig. 5f, the level lines are weighted with the average of gradient's magnitude. It thus provides a hierachical segmentation represented by its saliency map.

In Fig. 6, we show an illustration of our method applied to an image of license plate. The simplified image reconstructed from the set of disjoint level lines and the saliency map are depicted respectively in Fig. 6b and in Fig. 6c. An optimal thresholding of the saliency map is illustrated in Fig. 6d.

CONCLUSION AND PERSPECTIVES

In this paper, we have presented an efficient algorithm that selects disjoint level lines from a tree of shapes. The selection process is based on the incremental computation of the smallest enclosing shapes. The method extracts a set of disjoint level lines from meaningful ones to meaningless ones, that briefly reflect the main structure of the tree of shapes. Not only does it allow to analyse the tree structure, but it can also be a start for more advanced tree-based processing methods. Indeed, we have shown that it provides a relevant simplification of the images where the meaningful objects are separated one from each other and also provides a saliency map without transition problem around image contours. However, in natural images, the meaningful level lines are not all strictly separated, so in the future, we would like to relax this condition and tolerate a small superposition between the selected level lines. As a second future work, we aim at applying level lines shortening that yields an image curvature microscope [START_REF] Ciomaga | Level lines shortening yields an image curvature microscope[END_REF] before extracting disjoint level lines. Another major perspective is to extend this method to 3D images and the tree of shapes for color images [START_REF] Carlinet | A morphological tree of shapes for color images[END_REF][START_REF] Carlinet | Getting a morphological tree of shapes for multivariate images: paths, traps, and pitfalls[END_REF].

Fig. 2 :

 2 Fig. 2: An example of the tree-based image representations relying on threshold decomposition.

Fig. 3 :

 3 Fig.3: An example of Khalimsky's grid by materializing the points in a 2D image with 0-faces (blue disks), 1-faces (green strips), and 2-faces (red squares).

1 7 q 1 :

 71 COMPUTE SES(parent, S, depth) 2 foreach x in K Ω do SES(x) ← getCanonical(x) 3 foreach 2-face x in reverse order of S do 4 foreach 0 and 1-face e in x do 5 if depth(e) < depth(SES(x)) then 6 SES(x) ← getCanonical(e) ← parent(x) 8 if depth(SES(x)) < depth(SES(q)) then 9 SES(q) ← SES(x) 10 return SES Algorithm Computation of SES.

1 6 y 8 y

 68 SELECT LEVEL LINES(parent, SES, Order) 2 foreach x in K Ω do status(x) ← Null 3 S ′ = ∅ 4 foreach canonical element x in Order do 5 if status(x) = Unactive then ← parent(x) 7 while y = SES(x) and status(y) = Active do ← parent(y) 9 if y = SES(x) then 10 status(x) ← Active 11 S ′ ← S ′ ∪ {x} 12 y ← parent(x) 13 while y = SES(x) do 14 status(y) ← Unactive 15 y ← parent(y) 16 return S ′ Algorithm 2: Disjoint level lines selection algorithm.

Fig. 4 :

 4 Fig. 4: An illustration of the disjoint level lines selection with different orders. From top to down: original image, bottomup, top-down, average of gradient's magnitude ordering. Left: grayscale image; Right: corresponding randomly colorized image.

  (a) Input image f . (b) Simplified image f ′ . (c) Random color of f . (d) Random color of f ′ . (e) Segmented objects. (f) Saliency map.

Fig. 5 :

 5 Fig. 5: An illustration of the disjoint level lines selection applied to an image of bee board.

Fig. 6 :

 6 Fig. 6: An illustration of the disjoint level lines selection applied to an image of license plate.


Roland Levillain

Thierry Géraud
email: thierry.geraud@lrde.epita.fr

Laurent Najman
email: laurent.najman@esiee.fr

Edwin Carlinet

Practical Genericity: Writing Image Processing Algorithms Both Reusable and Efficient

Keywords: Image Processing, Software Reusability, Computational Efficiency, Generic Programming, Reproducible Research

An important topic for the image processing and pattern recognition community is the construction of open source and efficient libraries. An increasing number of software frameworks are said to be generic: they allow users to write reusable algorithms compatible with many input image types. However, this design choice is often made at the expense of performance. We present an approach to preserve efficiency in a generic image processing framework, by leveraging data types features. Variants of generic algorithms taking advantage of image types properties can be defined, offering an adjustable trade-off between genericity and efficiency. Our experiments show that these generic optimizations can match dedicated code in terms of execution times, and even sometimes perform better than routines optimized by hand.

Introduction

Like many numerical computing fields of computer science, Image Processing (IP) faces two kinds of issues which are hard to solve at the same time. On the one hand, the spectrum of data to process is very wide: IP methods can be used with "classical" 2-dimensional images with square pixels arranged along a regular grid, containing binary, gray-level, color, vectorial or even tensorial values; the dimension may also vary: 1D (signals), 3D (volumes), 2D+t (sequences of images) 3D+t (sequences of volumes); the domain of the image itself may be non-regular: images can be defined on mathematical structures such as graphs, simplicial complexes or topological maps.

The ability to handle so many data structures depends on the genericity of the theoretical framework and on the corresponding software tools as far as implementation is concerned. A generic algorithm is an algorithm that can be applied to a variety of inputs [START_REF] Géraud | Obtaining genericity for image processing and pattern recognition algorithms[END_REF], as opposed to a specific algorithm, which applies to a unique data type [START_REF] Dobie | Data structures for image processing in C[END_REF]. Generic Programming (GP) is a field of computer (c) On a 3D surface mesh (implemented as a simplicial complex).

Fig. 1: Results of the same morphological segmentation processing chain code applied to different input image types (a 2D square-pixel image, a graph, and a 3D surface); this particular example of a generic processing chain is detailed in [START_REF] Levillain | Writing reusable digital topology algorithms in a generic image processing framework[END_REF]. science (and a programming paradigm) which is concerned with issues of genericity in software 1 . Several software projects dedicated to IP rely on GP: Vigra [START_REF] Köthe | Reusable software in computer vision[END_REF][START_REF] Köthe | STL-style generic programming with images[END_REF], ITK [START_REF] Ibáñez | The ITK Software Guide[END_REF], Morph-M [START_REF]Centre for Mathematical Morphology: Morph-M: Image processing software specialized in mathematical morphology[END_REF], GIL [START_REF]Adobe: Generic Image Library (GIL)[END_REF]. Other noteworthy applications of GP to scientific computing include computational geometry [START_REF] Fabri | On the design of CGAL a computational geometry algorithms library[END_REF] and graph theory [START_REF] Siek | The Boost Graph Library[END_REF].

One of the incentives behind GP is the reusability of software, e.g. minimizing the cost of using existing algorithms with new data structures and vice versa. Such a cost may lie in bad software engineering practice such as duplicating code to process various types of inputs; or weakening typing (and most likely impact run-time efficiency) through generalization (e.g. using double as input type to also allow the use of bools, ints, or other types of values). Even well-engineered solutions, such as object-oriented ones, may harbour costs: dynamic dispatch (virtual methods) can become a prohibitively expensive feature w.r.t. execution times [START_REF] Géraud | Applying generic programming to image processing[END_REF].

On the other hand, many IP issues and applications involve large data sets (either numerous or voluminous) or make use of complex techniques requiring intensive computations. In both cases, practical software solutions have to meet efficiency constraints, especially regarding processing speed. Unfortunately, genericity and efficiency are often conflicting: efficient IP software is mostly dedicated to specific image types, methods or domains (and is therefore not generic). Conversely, most generic frameworks do not provide performances as good as specific ones.

In this paper, we investigate the issue of genericity versus efficiency with respect to IP, from the algorithmic point of view (i.e., we do not cover hardwarebased or assembly-level optimizations). We first explain how GP can be applied to IP, and what are the benefits of this approach (Section 2). Performance considerations in a generic context are addressed in Section 3. We study the causes of the opposition between genericity and efficiency and propose a trade-off as an answer to this problem, generic optimizations. The idea of generic optimization is carried further in Section 4, by not only acting on algorithm implementations, but by also working on input data types, so as to significantly increase performances at the expense of some genericity. Results of numerical experiments are presented and discussed in Section 5. The proposal of this paper is illustrated with an example from the field of Mathematical Morphology. The underlying ideas, are however very general, independent from the context (platform, compiler, etc.), and applicable to virtually any IP algorithm.

Genericity in Image Processing

A generic IP framework provides algorithms and data structures that have a single implementation and that can be used together to virtually implement any (valid) combination. This approach avoids duplication of similar code and enables a true reusability of algorithms on any compatible data structure (e.g., image type) without suffering from the combinatorial explosion. For example, the result images (right-hand column) presented in Figure 1 have been obtained using the same segmentation code on three images of very different types.

The simplest (though very limited) form of genericity in IP consists in parameterizing the type of values contained in images [START_REF] Ritter | Image Algebra: an overview[END_REF][START_REF] Kohl | The development of the Image Understanding Environment[END_REF], similarly to C++ standard containers [START_REF]Programming languages -C++[END_REF]. However, genuine genericity is based on the GP paradigm.

The key idea of this approach is to design the target framework using concepts, representing abstract entities of the domain (here, image processing) [START_REF] Levillain | Why and how to design a generic and efficient image processing framework: The case of the Milena library[END_REF]. Concepts define relationships between the corresponding entity (e.g. an image type) and other elements (e.g. image point type, image value type), as well as the minimal set of provided services (e.g. obtaining the value associated to a point). Generic algorithms are then written using concepts instead of specific concrete data types. As they expose no detail on the manipulated data types, such algorithms are generic implementations not tied to a specific input type. We have successfully applied this approach to Mathematical Morphology (MM) [START_REF] Darbon | Generic implementation of morphological image operators[END_REF][START_REF] Levillain | Milena: Write generic morphological algorithms once, run on many kinds of images[END_REF], discrete geometry [START_REF] Levillain | Writing reusable digital topology algorithms in a generic image processing framework[END_REF], as well as other fields of IP. In particular, one of the image dilation(const image& input) { image output(input.nrows(), input.ncols()); for (unsigned r = 0; r < input.nrows(); ++r) for (unsigned c = 0; c < input.ncols(); ++c) { unsigned char sup = input(r,c); for_all(q) if (input.has(q)) sup.take(input(q)); output(p) = sup.to_result(); } return output; } Algorithm 1.2: Generic dilation implementation [START_REF] Levillain | Milena: Write generic morphological algorithms once, run on many kinds of images[END_REF]. benefits of the generic approach is to enable user to try and experiment with uncommon and new data structures such as the Tree of Shapes [START_REF] Géraud | A quasi-linear algorithm to compute the tree of shapes of n-D images[END_REF][START_REF] Crozet | A first parallel algorithm to compute the morphological tree of shapes of nD images[END_REF][START_REF] Xu | Context-based energy estimator: Application to object segmentation on the tree of shapes[END_REF][START_REF] Xu | Salient level lines selection using the Mumford-Shah functional[END_REF], easily and rapidly.

if (r != 0 && input(r-1,c) > sup) sup = input(r-1,c); if (r != input.nrows()-1 && input(r+1,c) > sup) sup = input(r+1,c); if (c != 0 && input(r,c-1) > sup) sup = input(r,c-1); if (c != input.ncols()-1 && input(r,c+1) > sup) sup = input(r,
Let us illustrate the topic using a simple MM algorithm: a morphological dilation using a flat structuring element [START_REF] Géraud | Mathematical Morphology-From Theory to Applications[END_REF]. Algorithm 1.1 shows a simple C++ implementation of this filter. It however includes several implementation details which bounds the routine to specific inputs (2D image on a regular grid, having scalar values compatible with unsigned char). Moreover the (4-connected) structuring element cannot be changed. Therefore, we cannot use this algorithm to process, e.g., a 3D image composed of RGB color values and using 6-connectivity.

On the other hand, Algorithm 1.2 proposes a generic version of the previous morphological dilation. Input (resp. image and structuring element) types are now parameters of the algorithm (resp. I, an image type, and W, a window type); loops on vertical and horizontal ranges have been replaced by a single object p traversing the domain of the image input, called a site iterator on input's domain; likewise, members of the (previously hard-coded) structuring element w.r.t. p are no longer explicitly mentioned, as a another iterator q on the window (a qiter ) is used for this purpose; and instead of a hand-made computation of the maximum value, an accumulator object is used to iteratively compute the supremum from the values within the sliding window.

Efficiency vs Genericity Trade-off

Following the ideas expressed the previous section, we have designed and implemented a generic C++ IP library, Milena3 , which is part of the Free Software IP platform Olena [START_REF] Lrde | The Olena image processing platform[END_REF]. In addition to providing generic algorithms and data structures, Milena offers an extensible mechanism to implement efficient and reusable variants of some routines. In this section, we show how to implement efficient algorithms displaying a generic nature and how to wholly integrate them in the generic framework so that their use can be made transparent.

The Cost of Abstraction

Figures from Table 1 exhibit an important run time overhead in the generic case (Algorithm 1.2), which is about ten times longer to execute than the non-generic one (Algorithm 1.1). This is not a consequence of the GP paradigm per se. It is actually because of the highly abstract style of Algorithm 1.2, which in return makes the routine very versatile with respect to the context of use. The nongeneric version is faster than the generic one because it takes advantage of known features of the input types. For instance the structuring element is "built in the function", whereas it is an object taken as a generic input in Algorithm 1.2. Therefore its contents and size are constant and known at compile-time. Such implementation traits convey useful static information that compilers can leverage to optimize code. Hence, what prevents a code from being generic seems to be the condition to generate fast code: implementation details.

Generic Optimizations

The trade-off between genericity and efficiency admittedly depends on the level of details, but these two aspects are not entirely antagonistic: by carefully choosing the amount of specific traits used in an algorithm, one can create intermediate variants showing good run-time performance while keeping many generic features.

For instance, a means to speed up Algorithm 1.2 is to avoid using site iterators to browse the domain common to the input and output images. In Milena, site iterators can be automatically converted into sites (points), that is, locations in the domain of one (or several) image(s). Such location information is not tied to a given image: in the case of a regular 2D image, a site point2d(row, column) template <typename I, typename W> I dilation(const I& input, const W& win) { I output; initialize(output, input); // Iterator on the pixels of 'input'. mln_pixter(const I) pi(input); // Iterator on the pixels of 'output'. mln_pixter(I) po(output); // Iterator on the neighbor pixels of 'pi'. mln_qixter(const I, W) q(pi, win); for_all_2(pi, po) { accu::supremum<mln_value(I)> sup; for_all(q) sup.take(q.val()); po.val() = sup.to_result(); } return output; } Algorithm 1.3: Partially generic optimized dilation. is compatible with every regular, 2D, integer coordinate-based domain of the library (including toric spaces, non-rectangular 2D subspaces of Z 2 , etc.). This is why iterator p is used to refer to the same location in both input and output in Algorithm 1.2.

The price to pay for such a general expression is usually a run-time overhead: computations have to be performed each time a site iterator is used to access data from an image. However, this flexibility is not always needed when the data to process exhibit certain properties. For instance, an image whose values are stored in a contiguous, linear memory space, can be browsed using a pointer, directly accessing values in a sequential manner using their memory addresses, instead of computing a location at each access. In Milena, such pointers are encapsulated in small objects called pixel iterators or pixters where a pixel refers to an image's (site, value) pair. A pixter is bound to one image and cannot be used to iterate any other image. Pixters can also be used to browse spatiallyinvariant structuring elements (windows) as long as the underlying image domain is regular.

Algorithm 1.3 shows a reimplementation of Algorithm 1.2 where site iterators have been replaced by pixel iterators. The code is very similar, except that images input and output are now browsed with two different pixel iterators, each of them holding a pointer to the corresponding data. Such an implementation of the morphological dilation is less generic than the one of Algorithm 1.2. Even so, it can still be used with a wide variety of image types, as long as their data present a regular organization, which comprises any-dimension classical images using a single linear buffer to store their values. Besides, it is compatible with any spatially-invariant structuring element (or in other words, any constant window). Thus it remains much more generic than Algorithm 1.1. As for efficiency, Algorithm 1.3 matches almost Algorithm 1.1 in terms of speed (see Table 1), so it is a good alternative to the generic dilation, when the trade-off between and efficiency can be shifted towards the latter.

The approach presented here can be applied to other algorithms of the IP literature for which optimized implementations have been proposed. These optimizations are in practice compatible with a large set of input types, so their implementations can be considered as generic optimizations since they are not tied to a specific type.

Extra Generic Optimizations

The approach exposed in this paper can be carried further to improve the efficiency of generic optimizations. The idea is to involve data structures in the optimization effort: instead of acting only on algorithms, we can implement new optimized variants by working on their inputs as well.

For instance, in place of a window containing a dynamic array of vectorsthe size of which is known at run -we can implement and use a static window containing a static array carrying the same data, but whose size and contents are known at compile time. Modern compilers make use of this additional information to perform efficient optimizations (e.g, replace the loop over the elements of the window by an equivalent unrolled code). In this particular case, the implementation requires only the creation of two new, simple data types (static window, static pixel iterator). No additional implementation of the dilation is needed: Algorithm 1.3 is already compatible with this new window type. The resulting code delivers run times which are not only faster than the non-generic version of Algorithm 1.1, but which may also be faster than a handmade, pointer-based optimized (hence even less generic) version of the dilation, as shown in the next section.

Results

Table 1 shows execution times of various implementations of the morphological dilation with a 4-connected structuring element (window) applied to images of growing sizes (512 × 512, 1024 × 1024 and 2048 × 2048 pixels) . Times shown correspond to 10 iterative invocations. Tests were conducted on a PC running Debian GNU/Linux, featuring an Intel Pentium 4 CPU running at 3.4 GHz with 2 GB of RAM clocked at 400 MHz, using the C++ compiler g++ (GCC) version 4.4.5 with optimization options '-03' and '-DNDEBUG'.

In addition to the implementations shown in this paper, an additional nongeneric version using pointer-based optimizations has been added to the test suite, so as to further compare non-generic code -mostly optimized by handand generic code -mostly optimized by the compiler.

The overhead of the most generic algorithm is important: about ten times longer than Algorithm 1.1. The highly adaptable code of Algorithm 1.2 is free of implementation detail that the compiler could use to generate fast code (image values access with no indirection, statically-known structuring element). Algorithm 1.3 proposes a trade-off between genericity and efficiency: it is about 30%

Conclusion

This paper proposes an approach to reconcile genericity and efficiency in IP software. The key idea relies on generic optimizations expressed as algorithm specializations of the general case for a subspace of the acceptable input types. The addition of less generic but more efficient versions of an algorithm should not alter the motivation for designing an IP framework as generic as possible. We believe the most generic version of an algorithm should always be defined first, and then complemented by faster implementations. Firstly, having a generic version of an algorithm means having (at least) one implementation for each valid input type. Secondly, generic implementations are usually simpler, shorter and faster to write, provided the framework features entities supporting a generic programming style. Finally, generic implementations constitute a good basis to implement specializations, as their codes often share a similar structure.

The results presented in this paper are representative of the general outcomes of our proposal and are essentially independent from the compiler or platform used. In addition, we have already applied this strategy and observed the same conclusions regarding many other (and also more complex) algorithms that the one shown in this paper. Finally, we are not aware of any similar initiative regarding the efficiency of algorithms in generic IP libraries.

The Milena library, used to implement this paper's examples, is available in the Olena platform, a Free Software project released under the GNU General Public License that can be downloaded from our Web site [START_REF] Lrde | The Olena image processing platform[END_REF], as part of a reproducible research effort [START_REF] Buckheit | WaveLab and reproducible research[END_REF][START_REF] Fomel | Guest editors' introduction: Reproducible research[END_REF][START_REF] Morel | The IPOL initiative: Publishing and testing algorithms on line for reproducible research in image processing[END_REF]. This library is also a proof of concept of

  (a) On a "classical" regular grid. (b) On an edged-valued graph.

  Algorithm 1.1: Non-generic dilation implementation.

	c+1);
	output(r, c) = sup;
	}
	return output;
	}
	template <typename I, typename W>
	I dilation(const I& input, const W& win) {
	I output; initialize(output, input);

// Iterator on sites of the domain of 'input'. mln_piter(I) p(input.domain()); // Iterator on the neighbors of 'p' w.r.t. 'win'. mln_qiter(W) q(win, p); for_all(p) { // Accumulator computing the supremum on 'win'. accu::supremum<mln_value(I)> sup;

Table 1 :

 1 Execution times of various dilation implementations. times slower than Algorithm 1.1, but is generic enough to work on many regular image types (as a matter of fact, the most common ones). The case of the dilation with a static window is even more interesting: reusing the same code (Algorithm 1.3) with a less generic input (a static window representing a fixed and spatially-invariant structuring element) makes the code twice faster, to the point that it outperforms the manually optimized pointer-based implementation. Therefore, having several implementations (namely Algorithms 1.2 and 1.3) is useful when flexibility and efficiency are sought.

	Implementation	Time (s) per image (px)
		512 2 1024 2	2048 2
	Non generic (Alg. 1.1)	0.10 0.39	1.53
	Non generic, pointer-based 2 0.07 0.33	1.27
	Generic (Alg. 1.2)	0.99 4.07	16.23
	Fast, partly generic (Alg. 1.3) 0.13 0.54	1.95
	Alg. 1.3 with a static window 0.06 0.28	1.03

In this paper, by "generic" we mean "generic with respect to programming", not "generic with respect to the approach used to solve an IP issue". See for instance the definition of GP by Jazayeri et al.[START_REF] Jazayeri | Report of the Dagstuhl seminar (98061) on generic programming[END_REF].

Our library is available online from http://olena.lrde.epita.fr

This 26-line implementation is not shown in this paper for space reasons.the work presented in this paper. It features a collection of different image types along with many generic and efficient algorithms.