Stavros D Nikolopoulos
email: stavros@cs.uoi.gr

Charis Papadopoulos
email: charis@cs.uoi.gr

On the number of spanning trees of K m n ± G graphs

Keywords: Kirchhoff matrix tree theorem, complement spanning tree matrix, spanning trees, Kn-complements, multigraphs

The Kn-complement of a graph G, denoted by Kn -G, is defined as the graph obtained from the complete graph Kn by removing a set of edges that span G; if G has n vertices, then Kn -G coincides with the complement G of the graph G. In this paper we extend the previous notion and derive determinant based formulas for the number of spanning trees of graphs of the form K m n ± G, where K m n is the complete multigraph on n vertices with exactly m edges joining every pair of vertices and G is a multigraph spanned by a set of edges of K m n ; the graph

n by adding (resp. removing) the edges of G. Moreover, we derive determinant based formulas for graphs that result from K m n by adding and removing edges of multigraphs spanned by sets of edges of the graph K m n . We also prove closed formulas for the number of spanning tree of graphs of the form K m n ± G, where G is (i) a complete multipartite graph, and (ii) a multi-star graph. Our results generalize previous results and extend the family of graphs admitting formulas for the number of their spanning trees.

Introduction

The number of spanning trees of a graph G, denoted by τ (G), is an important, well-studied quantity in graph theory, and appears in a number of applications. Most notable application fields are network reliability [START_REF] Colbourn | The combinatorics of network reliability[END_REF][START_REF] Myrvold | Uniformly-most reliable networks do not always exist[END_REF][START_REF] Petingi | On the characterization of graphs with maximum number of spanning trees[END_REF], enumerating certain chemical isomers [START_REF] Brown | Some methods for counting the spanning trees in labelled molecular graphs, examined in relation to certain fullerenes[END_REF], and counting the number of Eulerian circuits in a graph [START_REF] Harary | Graph theory[END_REF].

Thus, both for theoretical and for practical purposes, we are interested in deriving formulas for the number of spanning trees of a graph G, and also of the K n -complement of G; the K n -complement of a graph G, denoted by K n -G, is defined as the graph obtained from the complete graph K n by removing a set of edges (of the graph K n ) that span G; if G has n vertices, then K n -G coincides with the complement G of the graph G. Many cases have been examined depending on the choice of G. For example, there exist closed formulas for the cases where G is is a pairwise disjoint set of edges [START_REF] Weinberg | Number of trees in a graph[END_REF], a chain of edges [START_REF] Moon | Enumerating labeled trees[END_REF], a cycle [START_REF] Gilbert | Maximizing spanning trees in almost complete graphs[END_REF], a star [START_REF] O'neil | The number of trees in a certain network[END_REF], a multi-star [START_REF] Nikolopoulos | On the number of spanning trees of multi-star related graphs[END_REF][START_REF] Yan | A formula for the number of spanning trees of a multi-star related graph[END_REF], a multi-complete/star graph [START_REF] Chung | On the number of spanning trees of a multi-complete/star related graph[END_REF], a labelled molecular graph [START_REF] Brown | Some methods for counting the spanning trees in labelled molecular graphs, examined in relation to certain fullerenes[END_REF], and more recent results when G is a circulant graph [START_REF] Golin | Chebyshev polynomials and spanning tree formulas for circulant and related graphs[END_REF][START_REF] Zhang | The number of spanning trees in circulant graphs[END_REF], a quasithreshold graph [START_REF] Nikolopoulos | The number of spanning trees in Kn-complements of quasi-threshold graphs[END_REF], and so on (see Berge [START_REF] Berge | Graphs and hypergraphs[END_REF] for an exposition of the main results).

In this paper, we extend the previous notion and consider graphs that result from the complete multigraph K m n by removing multiple edges; we denote by K m n the complete multigraph on n vertices with exactly m edges joining every pair of vertices. Based on the properties of the Kirchhoff matrix, which permits the calculation of the number of spanning trees of any given graph, we derive a determinant based formula for the number of spanning trees of the graph K m n -G, where G is a subgraph of K m n , and, thus, it is a multigraph. Note that, if m = 1 then K m n -G coincides with the graph K n -G. We also consider graphs that result from the complete multigraph K m n by adding multiple edges. More precisely, we consider multigraphs of the form K m n + G that result from the complete multigraph K m n by adding a set of edges (of the graph K m n ) that span G. Again, based on the properties of the Kirchhoff matrix, we derive a determinant based formula for the number of spanning trees of the graph K m n + G. To the best of our knowledge, not as much seems to be known about the number τ (K m n + G). Bedrosian in [START_REF] Bedrosian | Generating formulas for the number of trees in a graph[END_REF] considered the number τ (K n + G) for some simple configurations of G, i.e., when G forms a cycle, a complete graph, or when its vertex set is quite small. More recently, Golin et al. in [START_REF] Golin | Chebyshev polynomials and spanning tree formulas for circulant and related graphs[END_REF] derive closed formula for the number τ (K n + G) using Chebyshev polynomials, introduced in [START_REF] Boesch | Spanning tree formulas and Chebyshev polynomials[END_REF], for the case where G forms a circulant graph.

We denote K m n ± G the family of graphs of the forms K m n + G and K m n -G, and derive a determinant based formula for the number τ (K m n ± G). Moreover, based on these results, we generalize our formulas and extend the family K m n ± G to the more general family of graphs F m n ± G, where F m n is the complete multigraph on n vertices with at least m ≥ 1 edges joining every pair of vertices.

Based on our results-that is, the determinant based formulas for the number of spanning trees of the family of graphs τ (K m n ± G), and using standard algebraic techniques, we generalize known closed formulas for the number of spanning trees of simple graphs of the form K n -G. In particular, we derive closed formulas for the number of spanning trees τ (K m n ± G), in the case where G forms (i) a complete multipartite graph, and (ii) a multi-star graph.

We point out that our proposed formulas express the number of spanning trees τ (K m n ±G) as a function of the determinant of a matrix that can be easily constructed from the adjacency relation of the graph G. Our results generalize previous results and extend the family of graphs of the form K m n ± G admitting formulas for the number of their spanning trees.

Preliminaries

We consider finite undirected simple graphs and multigraphs with no loops; the term multigraph is used when multiple edges are allowed in a graph. For a graph G, we denote by V (G) and E(G) the vertex set and edge set of G, respectively.

The multiplicity of a vertex-pair (v, u) of a graph G, denoted by ℓ G (vu), is the number of edges joining the vertices v and u in G. The minimum multiplicity among all the vertex-pairs of G is denoted λ(G) while Λ(G) is the largest such number. Thus, if λ(G) > 0, then every pair of vertices in G is connected with at least λ(G) edges; if Λ(G) = 1, then G contains no multiple edges, that is, G is a simple graph (note that a simple graph or a multigraph contains no loops). The degree of a vertex v of a graph G, denoted by d G (v), is the number of edges incident with v in G. The minimum degree among the vertices of G is denoted δ(G) while ∆(G) is the largest such number.

It is worth noting that a multigraph can also be viewed as a simple graph whose edges are labeled with

On the number of spanning trees of K m n ± G graphs 237 non-zero integers (weights) which correspond to the multiplicity of the edges. Thus, a multigraph can be thought as a weighted graph. In this paper we adopt the standard approach.

Let F be the family of complete multigraphs on n vertices with multiplicity at least m ≥ 1, and let F m n ∈ F. Then, λ(F m n ) = m, since F m n has at least m ≥ 1 edges joining every pair of its vertices. A complete multigraph on n vertices with exactly m edges joining every pair of its vertices is called m-complete multigraph and denoted by K m n . Thus, for the m-complete multigraph K m n we have that

λ(K m n ) = Λ(K m n ) = m and δ(K m n ) = ∆(K m n ) = (n -1)m. Note that,
V (G + H) = V (G) ∪ V (H) and vu ∈ E(G + H) ⇐⇒ vu ∈ E(G) or vu ∈ E(H).
By definition, both graphs G and H are subgraphs of

G + H. Moreover, if v, u ∈ V (G) ∩ V (H), then ℓ G+H (vu) = ℓ G (vu) + ℓ H (vu).
Let G and H be two multigraphs such that E(H) ⊆ E(G). The graph G -H is defined as the graph obtained from G by removing the edges of H.

Having defined the graphs G + H and G -H, it is easy to see that

F m n = K m n + H(F m n ) and K m n = F m n -H(F m n ). In general, H(F m n ) = F m n -K m n ; the equality holds if H(F m n ) has n vertices. The adjacency matrix of a multigraph G on n vertices, denoted by A(G), is an n × n matrix with diagonal elements A(G)[i, i] = 0 and off-diagonal elements A(G)[i, j] = ℓ G (v i v j ). The degree matrix of the multigraph G, denoted by D(G), is an n × n matrix with diagonal elements D(G)[i, i] = d G (v i )
and off-diagonal elements D(G)[i, j] = 0. Throughout the paper empty entries in matrices represent 0s.

For an n × n matrix M , the (n -1)-st order minor µ i j is the determinant of the (n -1) × (n -1) matrix obtained from M after having deleted row i and column j; the i-th cofactor equals µ i i . The Kirchhoff matrix L(G) (also known as the Laplacian matrix) for a multigraph G on n vertices is an n × n matrix with elements

L(G)[i, j] = d G (v i ) if i = j, -ℓ G (v i v j ) otherwise where d G (v i ) is the degree of vertex v i in the graph G and ℓ G (v i v j )
is the number of edges joining the vertices v i and v j in G. The matrix L(G) is symmetric, has nonnegative real eigenvalues and its determinant is equal to zero. Note that,

L(G) = D(G) -A(G).
The Kirchhoff matrix tree theorem [START_REF] Biggs | Algebraic Graph Theory[END_REF] is one of the most famous results in graph theory. It provides a formula for the number of spanning trees of a graph G in terms of the cofactors of G's Kirchhoff matrix; it is stated as follows:

Theorem 2.1 (Kirchhoff Matrix Tree Theorem [START_REF] Biggs | Algebraic Graph Theory[END_REF]) For any multigraph G with L(G) defined as above, the cofactors of L(G) have the same value, and this value equals the number of spanning trees τ (G) of the multigraph G.

The Kirchhoff matrix tree theorem provides a powerful tool for computing the number τ (G) of spanning trees of a graph G. For this computation, we first form the Kirchhoff matrix L(G) of the graph G and obtain the (n -1) × (n -1) matrix L i (G) from L(G) by removing its i-th row and column (arbitrary), and then compute the determinant of the matrix L i (G). We note that the determinant of L i (G) simply counts the number of spanning trees rooted at vertex v i ∈ V (G), which justifies the reason of removing the i-th row and the i-th column from L(G). A combinatorial proof of the Kirchhoff matrix tree theorem can be found in [START_REF] Chaiken | A combinatorial proof of the all minors matrix tree theorem[END_REF].

The number of spanning trees of a graph G can be computed directly (without removing any row or column) in terms of a matrix L ′ (G) similar to the Kirchhoff matrix L(G), which is associated with the graph G [START_REF] Temperley | On the mutual cancellation of cluster integrals in Mayer's fugacity series[END_REF], or, alternatively, it can be computed by defining a characteristic polynomial det (L(G) -xI) on L(G); the latter approach takes into account the computation of the eigenvalues of the matrix L(G) (see [START_REF] Boesch | Spanning tree formulas and Chebyshev polynomials[END_REF][START_REF] Brown | Some methods for counting the spanning trees in labelled molecular graphs, examined in relation to certain fullerenes[END_REF][START_REF] Cvetković | Spectra of graphs[END_REF][START_REF] Golin | Chebyshev polynomials and spanning tree formulas for circulant and related graphs[END_REF][START_REF] Petingi | On the characterization of graphs with maximum number of spanning trees[END_REF][START_REF] Zhang | The number of spanning trees in circulant graphs[END_REF]).

In our work, we express the number of spanning trees of a graph of the form K m n ± G, where G is a subgraph of K m n on p vertices, in terms of a p × p matrix B(G) associated with the graph G, and not in terms of an n × n matrix L(K m n ± G) associated with the whole graph K m n ± G.

3 The K m n ± G graphs

In this section, we consider graphs that result from the m-complete multigraph K m n by removing or/and adding multiple edges. We are interested in deriving determinant based formulas for the number of spanning trees of the graphs K m n -G and K m n + G, where G is a multigraph spanned by a set of edges S ⊆ E(K m n ). To this end, we define a parameter α as follows: (i) α = 1, for the case K m n + G, and (ii) α = -1, for the case K m n -G. In other words, α = ±1, according to K m n ± G. Based on the value of α we conclude with the following result.

Let G be a multigraph spanned by a set of edges of the graph K m n . We derive formulas for the number of spanning trees of the graph K m n ± G; the graph G has p ≤ n vertices and Λ(G) ≤ m. In order to compute the number τ (K m n -G) we will make use of Theorem 2.1. Thus, we consider the n × n Kirchhoff matrix L = L(K m n ± G), which has the form:

L = 2
m(n -1)• • • -m -m -m • • • -m . . . . . . . . . . . . . . . . . . . . . -m • • •m(n -1) -m -m • • • -m -m • • • -m m(n -1) + α • dG(v1) -m -α • ℓG(vj vi) -m • • • -m m(n -1) + α • dG(v2) . . . • • • . . . . . . -m • • • -m -m -α • ℓG(vivj ) m(n -1) + α • dG(vp) 3
where

d G (v i ) is the degree of the vertex v i ∈ G and ℓ G (v i v j )
is the multiplicity of the vertices v i and v j in G. The entries of the off-diagonal positions (n -p + i, n -p + j) of the matrix L are equal to

On the number of spanning trees of K m n ± G graphs 239

-m -α • ℓ G (v i v j ), 1 ≤ i, j ≤ p.
Note that, the first n -p rows and columns of L correspond to the n -p vertices of the set V (K m n ) -V (G) and, thus, they have degree m(n -1) in K m n ± G. Let L 1 be the (n -1) × (n -1) matrix obtained from L by removing its first row and column. Then, from Theorem 2.1 we have that

τ (K m n -G) = det (L 1 ).
In order to compute the determinant of the matrix L 1 , we add one row and one column to the matrix L 1 ; the resulting n × n matrix L ′ 1 has 1 in position (1, 1), -m in positions (1, j), 2 ≤ j ≤ n, and 0 in positions (i, 1), 2 ≤ i ≤ n. It is easy to see that, det(L ′ 1 ) = det(L 1 ). Thus, the n × n matrix L ′ 1 has the following form:

L ′ 1 = 2
1 -m • • • -m -m -m • • • -m 0 m(n -1) • • • -m -m -m • • • -m . . . . . . . . . . . . . . . . . . . . . . . . 0 -m • • • m(n -1) -m -m • • • -m 0 -m • • • -m m(n -1) + α • dG(v1) -m -α • ℓG(vj vi) 0 -m • • • -m m(n -1) + α • dG(v2) . . . . . . • • • . . . . . . 0 -m • • • -m -m -α • ℓG(vivj ) m(n -1) + α • dG(vp)
.

We denote by L ′′ 1 the matrix obtained from L ′ 1 after multiplying the first row of L ′ 1 by -1 and adding it to the next n -1 rows. Thus, the determinant of L ′′ 1 is equal to the determinant of L ′ 1 . Moreover it becomes:

det (L ′′ 1 ) = 1 -m • • • -m -m -m • • • -m -1 mn . . . . . . -1 mn -1 mn + α • d G (v 1 ) -α • ℓ G (v j v i ) -1 mn + α • d G (v 2 ) . . . . . . -1 -α • ℓ G (v i v j ) mn + α • d G (v p )
,

where the entries of the off-diagonal positions (n -p + i, n -p + j) of the matrix

L ′′ 1 are equal to -α • ℓ G (v i v j ), 1 ≤ i, j ≤ p.
Note that, the first n -p rows of the matrix L ′′ 1 have non-zero elements in positions (1, i) and (i, i), 2 ≤ i ≤ n -p. We observe that the sum of all the elements on each row of L ′′ 1 , except of the first row, is equal to mn -

1; recall that, d G (v i ) = 1≤j≤p ℓ G (v i v j ), for every v i ∈ V (G).
Thus, we multiply each column of matrix L ′′ 1 by 1 mn and add it to the first column, and we obtain:

det (L ′′ 1 ) = 1 n -m • • • -m -m -m • • • -m 0 mn . . . . . . 0 mn 0 mn + α • d G (v 1 ) -α • ℓ G (v j v i ) 0 mn + α • d G (v 2 ) . . . . . . 0 -α • ℓ G (v i v j ) mn + α • d G (v p ) (2)
= m • (mn) n-p-2 • det (B),
where B = mnI p + αL(G) is a p × p matrix; recall that, L(G) is the Kirchhoff matrix of the multigraph G and thus, L(G) = D(G) ± A(G), where D(G) and A(G) are the degree matrix and the adjacency matrix of G, respectively. Concluding, we obtain the following result.

Theorem 3.1 Let K m n be the m-complete multigraph on n vertices, and let G be a multigraph on p vertices such that

V (G) ⊆ V (K m n ) and E(G) ⊆ E(K m n ). Then, τ (K m n ± G) = m • (mn) n-p-2 det(B), where B = mnI p + α • L(G) is a p × p matrix, α = ±1 according to K m n ± G, and L(G) is the Kirchhoff matrix of G.
We note that, for simple graphs K n and G in case K n -G, Theorem 3.1 has been stated first by Moon in [START_REF] Moon | Enumerating labeled trees[END_REF] and numerous authors in various guises used it as a constructive tool to obtain formulas for the number of spanning trees of graphs of the type K n -G.

In the previous theorem the graph G is a subgraph of K m n , and, thus, it has multiplicity Λ(G) ≤ m. It follows that the graph K m n + G has multiplicity Λ(K m n + G) ≤ 2m. However, we can relax the previous restriction in the case of the graph K m n + G. It is not difficult to see that for the matrix L of Equation (1) we have λ(G) ≥ 0. Thus, we can define the graph G to be a multigraph on p vertices such that V (G) ⊆ V (K m n ). The following theorem holds. Theorem 3.2 Let K m n be the m-complete multigraph on n vertices, and let G be a multigraph on p vertices such that

V (G) ⊆ V (K m n ). Then, τ (K m n + G) = m • (mn) n-p-2 • det(B)
, where B = mnI p + L(G) is a p × p matrix, and L(G) is the Kirchhoff matrix of G.

The F m n ± G graphs

In this section we derive determinant based formulas for the number τ (F m n ± G), where F m n is a complete multigraph and G is a subgraph of F m n . We first take into consideration the graph K m n + G 1 -G 2 and derive a determinant based formula for the number τ (K m n + G 1 -G 2 ), and, then, we derive a formula for the number τ (F m n ± G) using the graph

K m n + G 1 -G 2 and a characteristic graph H(F m n ).
On the number of spanning trees of K m n ± G graphs 241

4.1 The case K m n + G 1 -G 2
Here, we consider graphs that result from the m-complete multigraph K m n by adding multiple edges of a graph G 1 and removing multiple edges from a graph G 2 . Let G 1 be a multigraph on p 1 vertices, such that V (G 1 ) ⊆ V (K m n ), and let G 2 be a multigraph on p 2 vertices, such that

V (G 2 ) ⊆ V (K m n ) and E(G 2 ) ⊆ E(K m n + G 1 )
. We next focus on the graph K m n + G 1 -G 2 , which is obtained from the mcomplete multigraph K m n by first adding the edges of the graph G 1 and then removing from the resulting graph K m n + G 1 the edges of G 2 ; that is,

K m n + G 1 -G 2 = (K m n + G 1 ) -G 2 . It is worth noting that (K m n + G 1 ) -G 2 = K m n + (G 1 -G 2 ), since the graph G 2 is not necessarily a subgraph of G 1 . Moreover, V (G 1 ) = V (G 2 ).
For the pair of multigraphs (G 1 , G 2 ) we define the union-stable graphs G * 1 and G * 2 of (G 1 , G 2 ) as follows: G * 1 is the multigraph that results from G 1 by adding in V (G 1 ) the vertices of the set

V (G 2 ) - V (G 1 ) and G * 2 is the multigraph that results from G 2 by adding in V (G 2 ) the vertices of the set V (G 1 ) - V (G 2 ). Thus, by definition V (G * 1 ) = V (G * 2 ). By definition, both G * 1 and G * 2 are multigraphs on p = |V (G 1 ) ∪ V (G 2 )
| vertices, with at least p -p 1 and p -p 2 isolated vertices, respectively. Since

K m n + G * 1 -G * 2 = K m n + G 1 -G 2 , we focus on the graph K m n + G * 1 -G * 2 . Based on Theorem 2.1, we construct the n × n Kirchhoff matrix L = L(K m n + G * 1 -G * 2 ) associated with the graph K m n + G * 1 -G * 2 ;
it is similar to that of the case of K m n ± G. The difference here is the p × p submatrix which is formed by the last p rows and columns of L, where p = |V (G 1 ) ∪ V (G 2 )|. More precisely, the matrix L has the following form:

L =                  m(n -1) -m • • • -m -m -m • • • -m -m m(n -1) • • • -m -m -m • • • -m . . . . . . . . . . . . . . . . . . . . . . . . -m -m • • • m(n -1) -m -m • • • -m -m -m • • • -m B ′ (G * 1 , G * 2 )[1, 1] B ′ (G * 1 , G * 2 )[j, i] -m -m • • • -m B ′ (G * 1 , G * 2 )[2, 2] . . . . . . • • • . . . . . . -m -m • • • -m B ′ (G * 1 , G * 2 )[i, j] B ′ (G * 1 , G * 2 )[p, p]                  where the p × p submatrix B ′ (G * 1 , G * 2 ) has elements B ′ (G * 1 , G * 2 )[i, j] = m(n -1) + d G * 1 (v i ) -d G * 2 (v i ) if i = j, -m -ℓ G * 1 (v i v j ) + ℓ G * 2 (v i v j ) otherwise. We note that ℓ G * 1 (v i v j ) and ℓ G * 2 (v i v j ) are the number of edges of the vertices v i and v j in G * 1 and G * 2 , respectively. The entries d G * 1 (v i ) and d G * 2 (v i ), 1 ≤ i ≤ p, are the degrees of vertex v i of G * 1 and G * 2 , respectively. Note that, V (G * 1 ) = V (G * 2 )
. It is straightforward to apply a technique similar to that we have applied for the computation of the determinant of the matrix L ′ 1 in the case of K m n -G. Thus, in the case of

K m n + G * 1 -G * 2 , the determinant resulting (p + 1) × (p + 1) matrix B ′ (G) has 1 in position (1, 1), α in positions (1, j), 2 ≤ j ≤ p + 1,
and 0 in positions (i, 1), 2 ≤ i ≤ p + 1; recall that, α = ±1. Thus, the resulting matrix B ′ (G) has the following form:

B ′ (G) =         1 α α • • • α M 1 -α • • • -α -α M 2 • • • -α . . . . . . . . . . . . -α -α • • • M k         , where the diagonal m i × m i submatrices M i have diagonal elements mn + α • (p -m i ), 1 ≤ i ≤ k. Note that, det (B(G)) = det (B ′ (G)).
In order to compute the determinant of the matrix B ′ (G) we add the first row to the next p rows. Let B ′′ (G) be the resulting matrix. It follows that the determinant of B ′ (G) is equal to the determinant of B ′′ (G). We multiply the 2, 3, . . . , p + 1 columns of the matrix B ′′ (G) by -1/(mn + α • p) and add them to the first column; note that, the sum of each row of the matrix B ′′ (G) is equal to mn + α • p. Thus, the determinant of matrix B ′′ (G) becomes:

det (B ′′ (G)) = 1 -α•p mn+α•p α α • • • α M ′ 1 M ′ 2 . . . M ′ k = mn mn + α • p • det (M ′ 1 ) • det (M ′ 2 ) • • • • • det (M ′ k ), (3)
where the m i × m i submatrices M ′ i , 1 ≤ i ≤ k, have the following form:

M ′ i = mn + α • (p -m i ) + α α • • • α α mn + α • (p -m i ) + α • • • α . . . . . . . . . . . . α α • • • mn + α • (p -m i ) + α
For the determinant of matrix M ′ i we multiply the first row by -1 and add it to the next m i -1 rows. Then, we add the columns of matrix M ′ i to the first column. Observing that mn

+ α • (p -m i ) + α • m i = mn + α • p, we obtain det (M ′ i ) = (mn + α • p) • (mn + α • (p -m i )) mi-1
. Thus, from Equation (3) we have the following result.

Theorem 5.1 Let G = K m1,m2,...,m k be a complete multipartite graph on p = m 1 + m 2 + • • • + m k vertices. Then, τ (K m n ± G) = m • (mn) n-p-1 (mn + α • p) k-1 k i=1 (mn + α • (p -m i )) mi-1 ,
On the number of spanning trees of K m n ± G graphs 245

where p ≤ n and α = ±1 according to K m n ± G.

Remark 5.1. The class of complete multipartite graphs contains the class of c-split graphs (complete split graphs); a graph is defined to be a c-split graph if there is a partition of its vertex set into a stable set S and a complete set K and every vertex in S is adjacent to all the vertices in K [START_REF] Golumbic | Algorithmic graph theory and perfect graphs[END_REF]. Thus, a c-split graph G on p vertices and V (G) = K +S is a complete multipartite graph K m1,m2,...,m k with m 1 = |S|, m 2 = m 3 = • • • = m k = 1 and k = |K| + 1. A closed formula for the number of spanning trees of the graph K n -G was proposed in [START_REF] Nikolopoulos | The number of spanning trees in Kn-complements of quasi-threshold graphs[END_REF], where G is a c-split graph.

Let G be a c-split graph on p vertices and let V (G) = K + S. Then, from Theorem 5.1 we obtain that the number of spanning trees of the graphs K m n ± G is given by the following closed formula:

τ (K m n ± G) = m • (mn) n-p-1 (mn + α • |K|) |S|-1 (mn + α • p) |K| ,
where p = |K| + |S| and p ≤ n. ✷ It has been proved [START_REF] Chung | On the number of spanning trees of a multi-complete/star related graph[END_REF][START_REF] Nikolopoulos | On the number of spanning trees of multi-star related graphs[END_REF][START_REF] Yan | A formula for the number of spanning trees of a multi-star related graph[END_REF] that the number of spanning trees of the graph K n -G is given by the following closed formula:

Multi-star graphs

τ (K n -G) = n n-p-2 (n -1) p-r 1 + r i=1 1 q i -1 • r i=1 (q i -1) ,
where

q i = n -(r -1 + b i ) -bi n-1 .
In this section, based on Theorem 3.1, we generalize the previous result by deriving a closed formula for the number of spanning trees of the graphs K m n ± G, where G is a multi-star on p ≤ n vertices. Let K r be the complete graph of the multi-star graph G and let v 1 , v 2 , . . . , v r be its vertices. The vertex set consisting of the vertex v i and the b i vertices of degree one which are incident with vertex v i induces a star on b i + 1 vertices, 1 ≤ i ≤ r. We construct a (b i + 1) × (b i + 1) matrix M i which corresponds to the star with center vertex v i ; it has the following form:

M i =      mn + α -α mn + α -α . . . -α -α • • • mn + α • (r -1 + b i )     
, where α = ±1.

In order to compute the determinant of the matrix M i we first multiply the first row by -1 and add it to the next b i -1 rows. We then add the b i columns to the first column. Finally, we multiply the first column by α mn+α and add it to the last column. Thus, by observing that α 2 = 1, we obtain:

det (M i ) = (mn + α) bi • mn + α • (r -1 + b i ) - α 2 b i mn + α = (mn + α) bi • mn + α • (r -1 + b i ) - b i mn + α = (mn + α) bi • q i , where q i = mn + α • (r -1 + b i ) - b i mn + α , 1 ≤ i ≤ r. (4)
We are now in a position to compute the number of spanning trees τ (K m n ± G) using Theorem 3.1. Thus, we have

τ (K m n ± G) = m • (mn) n-p-2 • det (B(G)) (5)
where

B(G) =               M 1,1 -α M 2,2 -α . . . . . . M r,r -α -α mn + α • d G (v 1 ) -α • • • -α -α -α mn + α • d G (v 2 ) • • • -α . . . . . . . . . . . . . . . -α -α -α • • • mn + α • d G (v r )              
is a p × p matrix and M i,i is a submatrix which is obtained from M i by deleting its last row and its last column,

1 ≤ i ≤ r. The degrees of the vertex v i of K r is equal to d G (v i ) = r -1 + b i , 1 ≤ i ≤ r.
It now suffices to compute the determinant of the matrix B(G). Following a procedure similar to that we applied to the matrix M i , we obtain:

det (B(G)) = (mn + α) p-r • q 1 -α • • • -α -α q 2 • • • -α . . . . . . . . . . . .
-α -α • • • q r = (mn + α) p-r • det (D).

Recall that, q i = mn + α • (r -1 + b i ) -bi mn+α ; see Equation [START_REF] Boesch | Spanning tree formulas and Chebyshev polynomials[END_REF]. In order to compute the determinant of the r × r matrix D we first multiply the first row of D by -1 and add it to the r -1 rows. Then, we multiply column i by q1+α qi+α , 2 ≤ i ≤ r, and add it to the first column. Expanding in terms of the rows of matrix D, we have that

det (D) = 1 -α r i=1 1 q i + α • r i=1
(q i + α) .

Thus, substituting the value of det (D) into Equation ( 5), we obtain the following theorem.

On the number of spanning trees of K m n ± G graphs 247 (q i + α) , where p ≤ n, q i = mn + α • (r -1 + b i ) -bi mn+α and α = ±1 according to K m n ± G.

Concluding remarks

In this paper we derived determinant based formulas for the number of spanning trees of the family of graphs of the form K m n ± G, and also for the more general family of graphs F m n ± G, where K m n (resp. F m n ) is the complete multigraph on n vertices with exactly (resp. at least) m edges joining every pair of vertices and G is a multigraph spanned by a set of edges of K m n (resp. K m n ). Based on these determinant based formulas, we prove closed formulas for the number of spanning trees τ (K m n ± G), in the case where G is (i) a complete multipartite graph, and (ii) a multi-star graph.

In light of our results, it would be interesting to consider the problem of proving closed formulas for the number of spanning tree τ (K m n ± G) in the cases where G belongs to other classes of simple graphs or multigraphs. Moreover, instead of a multigraph G that we have taken into account, we do not know whether a generalazisation for an arbitrarily weighted graph G holds.

The problem of maximizing the number of spanning trees was solved for some families of graphs of the form K n -G, where G is a multi-star graph, a union of paths and cycles, etc. (see [START_REF] Chung | On the number of spanning trees of a multi-complete/star related graph[END_REF][START_REF] Gilbert | Maximizing spanning trees in almost complete graphs[END_REF][START_REF] Nikolopoulos | On the number of spanning trees of multi-star related graphs[END_REF][START_REF] Petingi | On the characterization of graphs with maximum number of spanning trees[END_REF]). Thus, an interesting open problem is that of maximizing the number of spanning trees of graphs of the form K m n ± G.

A

  multi-star graph, denoted by K r (b 1 , b 2 , . . . , b r ), consists of a complete graph K r with vertices labelled v 1 , v 2 , . . . , v r , and b i vertices of degree one, which are incident with vertex v i , 1 ≤ i ≤ r [7, 19, 25]; see Figure 1(b). Let G = K r (b 1 , b 2 , . . . , b r ) be a multi-star graph on p = r + b 1 + b 2 + • • • + b r vertices.

Theorem 5 . 2

 52 Let G = K r (b 1 , b 2 , . . . , b r ) be a multi-star graph on p = r + b 1 + b 2 + • • • + b r vertices. Then, τ (K m n ± G) = m • (mn) n-p-2 (mn + α) p-r 1 -α •

Acknowledgements

The authors would like to express their thanks to the anonymous referees whose suggestions helped improve the presentation of the paper.

of the matrix L ′ 1 of Equation [START_REF] Berge | Graphs and hypergraphs[END_REF] becomes

where the p × p submatrix B(G * 1 , G * 2 ) has elements

and obtain the following result. Theorem 4.1 Let K m n be the m-complete multigraph on n vertices, and let

,

and L(G * 2 ) are the Kirchhoff matrices of the union-stable graphs G * 1 and G * 2 of (G 1 , G 2 ), respectively.

The general case F m n ± G

Let F m n be a complete multigraph on n vertices and let G be a subgraph of F m n . We will show that the previous theorem provides the key idea for computing the number τ (F m n ± G), where F m n ± G is the graph that results from F m n by adding or removing the edges of G. Since λ(F m n ) > 0, we have that

The addition of the edges of G in the graph F m n , implies that

Thus, for the computation of the number τ (F m n + G) we can apply Theorem 3.2. On the other hand, in the case of removal the edges of G from the graph F m n , for the computation of the number τ (F m n -G) we can apply Theorem 4.1, since

Concluding we have the following result. where

Note that, we consider the graph F m n ± G and therefore G must be a subgraph of F m n . However in the case of the F m n + G graph, similar to Theorem 3.2, it is obvious that G can be a graph spanned by any set of edges joining the vertices of F m n .

Classes of graphs

In this section, we generalize known closed formulas for the number of spanning trees of families of graphs of the form K n -G. As already mentioned in the introduction there exist many cases for the τ (K n -G), depending on the choice of G. The purpose of this section is to prove closed formulas for τ (K m n ±G), by applying similar techniques to that of the case of K n -G. Thus we derive closed formulas for the number of spanning trees τ (K m n ±G), in the cases where G forms (i) a complete multipartite graph, and (ii) a multi-star graph.

Complete multipartite graphs

A graph is defined to be a complete multipartite (or complete k-partite) if there is a partition of its vertex set into k disjoint sets such that no two vertices of the same set are adjacent and every pair of vertices of different sets are adjacent. We denote a complete multipartite graph on p vertices by K m1,m2,...,m k , where p = m 1 + m 2 + • • •+ m k ; see Figure 1(a). We note that the number of spanning trees of a complete multipartite graph has been considered by several authors in the past [START_REF] Clark | On the enumeration of spanning trees of the complete multipartite graph[END_REF][START_REF] Lewis | The number of spanning trees of a complete multipartite graph[END_REF].

Let G = K m1,m2,...,m k be a complete multipartite graph on p vertices. In [START_REF] O'neil | Enumeration of spanning trees in certain graphs[END_REF] it has been proved that the number of spanning trees of K n -G is given by the following formula:

where p is the number of vertices of G.

In this section, we extend the previous result by deriving a closed formula for the number of their spanning trees for the graphs K m n ± G, where G a complete multipartite graph on p ≤ n vertices. From Theorem 3.1, we construct the p × p matrix B(G) and add one row and column to the matrix B(G); the


Mathieu Dubois
email: mathieu.dubois@limsi.fr

Emmanuelle Frenoux
email: emmanuelle.frenoux@limsi.fr

Philippe Tarroux
email: philippe.tarroux@limsi.fr

Using n-grams models for visual semantic place recognition

Keywords: Semantic place recognition, Hidden Markov models, n-grams

The aim of this paper is to present a new method for visual place recognition. Our system combines global image characterization and visual words, which allows to use efficient Bayesian filtering methods to integrate several images. More precisely, we extend the classical HMM model with techniques inspired by the field of Natural Language Processing. This paper presents our system and the Bayesian filtering algorithm. The performance of our system and the influence of the main parameters are evaluated on a standard database. The discussion highlights the interest of using such models and proposes improvements.

INTRODUCTION

Semantic mapping (see (Nüchter and Hertzberg, 2008)) is a relatively new field in robotics which aims to give the robot a high-level, human-compatible, understanding of its environment in order to ease the integration of robots in daily environments, notably homes or workplaces. Such environments are usually composed of discrete places which correspond to different functions. For instance a house is usually made of different rooms and corridors used to move between them. Such places are called semantic places because they are defined in high-level human concepts as opposed to traditional low-level landmarks used in robot mapping.

In this context, it's important for the robot to be able to recognize in which place or category of places it lies. Those tasks are called respectively instance recognition and categorization. Semantic place recognition is then an important component of semantic mapping. Moreover the semantic category of a place can be used to foster object detection and recognition (giving priors on objects identity, location and scale) and to provide qualitative localization.

Different types of sensors have been employed for semantic place recognition. The first works in this domain used range sensors to discriminate places based on geometrical information. However the spatial configuration of two places of the same category (e.g. two kitchens) can be very different. Therefore geometrical information may not be useful for categoriza-tion. Vision is the modality of choice for semantic place recognition because it gives access to rich, allothetic information. Although there are multimodal approaches, our work focuses on visual place recognition.

In this article we will further develop an analogy between semantic place recognition and language modelling. This analogy allows to design efficient temporal integration methods i.e. to take several images into account in order to reduce ambiguity. More precisely, we will extend the Hidden Markov Model (HMM) formalism with n-grams models. Those models have been extensively used in Natural Language Processing (NLP) and efficient estimation techniques have been proposed. This paper aims to assess the use of such models in semantic place recognition. The goal is to compare this temporal integration method to previously proposed models. In particular we will study the influence of the length of the n-gram model and estimation procedure on performance.

The article is structured as follows. Section 2 presents related work. Our model and its links with language modelling are described in section 3. Section 4 presents our experiments and the results. Finally we conclude in section 5.

RELATED WORK

Some authors (see (Vasudevan et al., 2007)) use an object-based approach. In this case they employ a standard algorithm for object localization and recognition. Places are described by the frequency of objects found in them combined with constraints on their position. However, object categorization is still a difficult task and the position of objects can greatly vary from one environment to another. Therefore those approaches have not been used on large databases.

The vast majority of research on place recognition use techniques developed for visual scene classification. We can distinguish methods using global features (see (Torralba et al., 2003)) and methods using descriptors computed around interest points (see (Ullah et al., 2008)). (Filliat, 2008) uses the Bag-of-Words (BoW) model: local features are first clustered into a so-called dictionary of visual words learned by mean of a vector quantization algorithm. An image is represented by the distribution of visual words found in it. The major advantage is that the learning space is discretized but all geometrical information is lost.

Generally speaking using a single image or a single type of information is not enough for place recognition tasks. Therefore a lot of research has been conducted to disambiguate perception. (Pronobis and Caputo, 2007) use a confidence criterion to iteratively compute several cues from the same image until confidence in the classification is sufficiently high (or no more cues are available).

Another method to reduce ambiguity is to use several images to mutually disambiguate perception. In (Pronobis et al., 2010), the authors use a simple spatio-temporal accumulation process to filter the decision of a discriminative confidence-based place recognition system (which uses only one image to recognize the place). One problem with this method is that the system needs to wait some time before giving a response. Also, special care must be taken to detect places boundaries and to adjust the size of the bins. (Torralba et al., 2003) use a HMM where each place is a hidden state and the feature vector stands for the observation. The drawback is that the input space is continuous and high-dimensional. The learning procedure is then computationally expensive. (Ranganathan, 2010) uses a technique called Bayesian online change-point detection. The main idea is to detect abrupt changes in the parameters of the input's statistics caused by moving from one place to another. The main advantage is that the robot is able to learn in an unsupervised way but relies on the hypothesis that the shape of the distribution is the same for every place.

Several works (see (Wu et al., 2009;Guillaume et al., 2011;Dubois et al., 2011)) have combined global image description and vector quantization. In this case, each image is described by a single visual word. The sequence of images is then translated into a sequence of words. Such techniques allow to draw a parallel between place recognition and language modelling. (Wu et al., 2009) propose to use a HMM with discretized signatures. Temporal integration is performed with Bayesian filtering (see section 3). (Dubois et al., 2011) propose to use an extended model called auto-regressive HMM to take into account the dependence between images.

In this paper we push this idea a step further. The next section presents our models and its relations to the standard HMM model.

PLACE RECOGNITION WITH n-GRAMS

Our model is similar to the one described in (Guillaume et al., 2011;Dubois et al., 2011). Each image is described by a unique feature vector which is mapped to a given visual word thanks to a vector quantization algorithm (see section 3.3). The main novelty lies in the use of High-Order Hidden Markov Model (see section 3.1) and techniques for visual word selection (see 3.4).

High-Order Hidden Markov model

In HMMs the relationship between x t , the robot's knowledge of the world at time t, and z t , its perception is represented by figure 1(a). In the case of place recognition, the state is a discrete random variable which represents the place the robot is in at time t. In this model, each place c i ∈ C is modelled by the continuous probability distribution p (z t |x t = c i ). This formalism allows to efficiently estimate the a posteriori probability bel(x t ) = P (x t |z 1:t ) by a recursive equation (see (Wu et al., 2009)) given the discrete place transition probability distribution P (x t |x t-1 ) which encodes the topology of the environment.

It is assumed that the current observation depends only on the current hidden state i.e. that the state is complete. However, there is a huge semantic gap between the human notion of a place and what can be extracted from an image. Several authors have proposed extensions of the classic HMM to take into account long-term dependencies between observations (see [START_REF] Berchtold | High-order extensions of the double chain markov model[END_REF]Lee and Lee, 2006)). In this paper we will call this model High-Order Hidden Markov Model (HOHMM). In this case, the current knowledge x t depends on the last ℓ states x t-ℓ:t-1 . Similarly the current observation z t depends on x t and the n previous observations z t-n:t-1 (see figure 1(b)). In this paper we restrict ourselves to the case ℓ = 1. Therefore the state transition matrix is unchanged.

The a posteriori distribution bel (x t ) is given by: bel

(x t ) = p (z t |z t-n:t-1 , x t ) ∑ c i ∈C P (x t |x t-1 ) bel (x t-1 )
(1) The place model is given by the distribution p (z t |z t-n:t-1 , x t = c i ). This probability distribution may be very difficult to learn because it is continuous.

HOHMM and Visual Words

In order to simplify learning of the place model (Guillaume et al., 2011;Dubois et al., 2011) have proposed to use global image characterization in combination with vector quantization algorithms to discretize them. In this case the variable z t is reduced to a discrete random variable with a finite number of values {1, . . . , K} where K is the number of words in the dictionary.

In this case, the model of place c i is given by the discrete probability distribution P (z t |z t-n:t-1 , x t = c i ). In NLP, such a model is known as a n + 1-gram model because it uses n + 1 words. [START_REF] Chen | An empirical study of smoothing techniques for language modeling[END_REF] have shown that the estimation of the model from empirical data is an important factor. One problem is that even with a large training set, some sequences of words will not be observed in training data for a given class and therefore they will be assigned a null probability in this class' model. If such a sequence is observed in the testing set then the a posteriori probability of this class will be clamped to 0 due to equation 1. To avoid this problem, it is necessary to take some probability mass from the observed sequences and distribute it to unobserved sequences. Those techniques are called smoothing or discounting. We refer the reader to (Manning and Schütze, 1999) for a unified presentation of smoothing techniques. We use the SRILM toolkit to learn the n-grams models.

Image Characterization and Vector Quantization

To characterize the images we use the GIST descriptors (see (Torralba et al., 2003)) which is an efficient global image characterization. The image is divided into 4 × 4 subwindows (we use only the luminance channel) and filtered using a bank of Gabor filters (we use 4 scales and 6 orientations). The energy of the filter is then averaged on each subwindow for each scale and orientation. Finally the output is projected on the first 80 principal components which explains more than 99% of the variance. Thus this descriptor captures the most significant spatial orientations at a given scale.

The vector quantization algorithm used in this paper is the Self-Organizing Map (SOM) (see (Kohonen, 1990)). In the current set-up the training of the SOM is performed off-line on a set of randomly chosen images made of 1 /3 of the COLD DB. The number of neurons on the map sets the number of words in the visual dictionary which is an important parameter of the system. We use square maps parametrized by their length S (therefore K = S 2 ). In this paper we will use S = 10 and S = 20. Those values were selected because it has been shown that small maps have a good performance on categorization tasks while larger maps perform well for instance recognition (see (Guillaume et al., 2011)). Because the training algorithm is stochastic, the results vary from one SOM to another. Therefore for each size S, the results are averaged for 5 SOMs.

Visual Words Selection

The sampling rate of most databases is several Hertz. In this case, image at time t + 1 is not very different from image at time t and there is a high probability that they are described by close vectors and therefore by the same visual word. While this is a desirable feature of image description and vector quantization, this may be a problem for our method because the probability of seeing the same visual word than before will be very high. Therefore it might be interesting to use only a subset of the images (and then the words) for learning.

In order to evaluate this phenomenon we have computed the average number of consecutive timesteps which are characterized by the same visual word for the training sequence used in section 4. Results are given in table 1.

We will test three different strategies for selecting visual words. The first one is simply to sub-sample the input image i.e. to select 1 image out of s (s is the sub-sampling rate). This strategy will be called "subsample". The second strategy is to replace every sequence of m identical prototypes by a unique instance of this word (m is the compression rate). We will call this strategy "compress". The last strategy is to use the word at time t only if it is different than the word at time t -1. We will call this strategy "unique". Those strategies are simple and can be implemented online on a real robot with limited computational power.

In the next section we will present the experiment we carried out to study the use of this model for se-

X t-1 X t Z t-1 Z t (a) X t-ℓ X t-1 X t . . . Z t-n Z t-2 Z t-1 Z t . . .

EXPERIMENTAL RESULTS

Experimental Design

We use the COLD database (see (Ullah et al., 2008)) a standard database to evaluate vision-based place recognition systems. It consists of sequences acquired by a human-driven robot in different laboratories across Europe under different illumination conditions (night, cloudy, sunny). In each laboratory, two paths were explored (standard and extended). Each path was followed at least 3 times under each illumination condition. All the experiments were carried out with the perspective images.

Protocols proposed by (Ullah et al., 2008) uses only a few hundreds images per place which is not enough to robustly estimate the transition probabilities. Therefore we designed a new experiment to evaluate the interest of our method. We use only images acquired in Saarbruecken part B because other parts of the database are known to contain errors (e.g. missing places or labellisation errors) or are not complete (e.g. only one path was followed). There are five classes (see table 2). Training is performed with sequences number 1 and 2 from all the three illumination conditions. Similarly, testing is performed with sequence 3 from all the illumination conditions.

Following (Wu et al., 2009) we define the transition matrix as P (x t |x t-1 ) = p e if x t = x t-1 ; the rest of the probability mass is shared uniformly among all other transitions. We use p e = 0.99.

In order to test the influence of the n-gram order we have varied n between 1 and 6. Similarly we have tested the Lidstone-Laplace (LL) smoothing with parameter δ = 1 and the Witten-Bell (WB) smoothing.

The training set was too small to use the Knesser-Nay smoothing. In our experiments we use interpolated models (Manning and Schütze, 1999). We have tested several values of the sub-sampling rate: s = 1 (which has no effect), s = 3 and s = 5. We use m = 3 for the "compress" strategy. The "unique" strategy don't need any parameter.

Setting n = 1 with Lidstone-Laplace smoothing gives the same temporal integration method than in (Wu et al., 2009) (note that we don't use the same signature). Setting n = 2 with Lidstone-Laplace smoothing and without interpolation gives a system similar to (Dubois et al., 2011).

Results

Results are presented on figure 2. It must be noted that on this instance recognition task, a larger SOM gives better results. This is expected from the literature (see (Guillaume et al., 2011)). The second observation that could be made is that the word selection methods generally increase the results by several percent. This can be seen by the difference between the bar for s = 1 and other bars of the same group. The "subsample" strategy with s = 3 is rather efficient, sometimes increasing performance by 6%. Setting s = 5 generally gives less important increase. Performance decreases with S = 20 and LL smoothing. However, this strategy leads to the best results on the task for S = 20 and WB smoothing. The "compress" strategy is usually efficient except for s = 20 and WB smoothing. The "unique" strategy is always among the best choices and it's results are less sensitive to the n-gram order. Generally speaking, the effect of those strategies increase with n. Results with WB smoothing are generally a little bit better than with LL in particular for large n.

It is clear from the figure that using n = 2, i.e. to take into account the dependence on the last image, is a clear improvement over n = 1, i.e. the classical HMM. However using n-grams with n > 2 has little impact on performance. It should be noted that when s = 1, the performance drops when n > 2. With word selection, the performance can be high with large n. This seems to confirm the intuition behind the word selection techniques.

CONCLUSION

We have presented a new model of temporal integration using HOHMM for semantic place recognition which models the dependence between observations. We have shown that taking this dependence into account can lead to interesting gains in performance. However, contrary to what we expected, using larger n don't improve performance. The smoothing technique seems to have minor effect. This may be caused by the fact that we use relatively small training sets compared to the field of NLP where those techniques have been developed. Those results must take into account the fact that recognition rates are already quite high on the task studied here.

We have shown that simple methods to select important words could improve the results. Our results suggest that large n could be interesting if combined with good word selection techniques.

Future works will focus on the vector quantization process to learn better words. More sophisticated word selection techniques may also be useful. Finally we could also look for more discriminative descriptors.

  (a) The classical HMM model. (b) The HOHMM model (we only show nodes that have an influence on x t and z t ).

  S = 20, Witten-Bell smoothing Figure2: Results on the instance recognition task. The vertical axis is the correct recognition rate (in %). The horizontal axis is the value of n. Upper-row: results for S = 10. Lower row: results for S = 20. Left column: results for Lidstone-Laplace smoothing. Right column: results for Witten-Bell smoothing.

Table 1 :

 1 Average number of consecutive images represented by the same visual word for different SOM size.

mantic place recognition.

Table 2 :

 2 Number of images for each category in the training and testing sets. There are 11,380 training and 5,192 testing images

		Office	Corridor	Printer area	Toilets	Kitchen
	Training	1,375	4,464	1,190	3,272	1,079
	Testing	606	1,964	532	1,513	577

ACKNOWLEDGEMENTS

We thanks Thiago Fraga and Alexandre Allauzen for fruitful discussion and help with the SRILM toolkit.

ceedings of the 34th annual meeting on Association for Computational Linguistics. Dubois, M., Guillaume, H., Tarroux, P., and Frenoux, E. (2011). Visual place recognition using bayesian filtering with markov chains. In Proceedings of the European Symposium on Artificial Neural Networks (ESANN 2011). Filliat, D. (2008). Interactive learning of visual topological navigation. In Proceedings of the 2008 IEEE International Conference on Intelligent Robots and Systems (IROS 2008). Guillaume, H., Dubois, M., Tarroux, P., and Frenoux, E. (2011). Temporal Bag-of-Words: A Generative Model for Visual Place Recognition using Temporal Integration. In Proceedings of the International Conference on Computer Vision Theory and Applications (VIS-APP 2011). Kohonen, T. (1990). The self-organizing map. In Proceedings of the IEEE, volume 78, pages 1464-1480. Lee, L.-M. and Lee, J.-C. (2006). A study on high-order hidden markov models and applications to speech recognition. In Ali, M. and Dapoigny, R., editors, Advances in Applied Artificial Intelligence, volume 4031 of Lecture Notes in Computer Science. Manning, C. and Schütze, H. (1999). Foundations of statistical natural language processing. MIT Press. Nüchter, A. and Hertzberg, J. (2008). Towards semantic maps for mobile robots. Robotics and Autonomous Systems, 56(11):915-926. Pronobis, A. and Caputo, B. (2007). Confidence-based cue integration for visual place recognition. In Proccedings of the 2007 IEEE/RSJ International Conference on Intelligent Robots and Systems. Pronobis, A., Mozos, O. M., Caputo, B., and Jenseflt, P. (2010). Multi-modal semantic place classification. The International Journal of Robotics Research, 29(2-3):298-320. Ranganathan, A. (2010). PLISS: Detecting and labeling places using online change-point detection. In Proceedings of the 2010 Robotics: Science and Systems Conference (RSS 2010). Torralba, A., Murphy, K. P., Freeman, W. T., and Rubin., M. A. (2003). Context-based vision system for place and object recognition. In Proceedings of the Nineth IEEE International Conference on Computer Vision (ICCV 2003), volume 1, pages 273-280. Ullah, M. M., Pronobis, A., Caputo, B., Luo, J., Jensfelt, P., and Christensen, H. I. (2008). Towards robust place recognition for robot localization. In Proceedings of the IEEE International Conference on Robotics and Automation (ICRA 2008), Pasadena, USA. Vasudevan, S., Gachter, S., Nguyen, V., and Siegwart, R. (2007). Cognitive maps for mobile robots-an object based approach. Robotics and Autonomous Systems, 55(5):359-371. Wu, J., Christensen, H., and Rehg, J. (2009). Visual place categorization: Problem, dataset, and algorithm. In IEEE/RSJ International Conference on Intelligent Robots and Systems, 2009 (IROS 2009).


Michel Vacher

Pedro Chahuara

Benjamin Lecouteux

Franc ¸ois Dan Istrate
email: dan.istrate@esigetel.fr

François Portet

Thierry Joubert
email: thierry.joubert@theoris.fr

Mohamed El

Amine Sehili

Brigitte Meillon

Nicolas Bonnefond

Sébastien Fabre
email: sebastien.fabre@technosens.fr

Thierry Portet

Mohamed Joubert

Brigitte Sehili

Nicolas Meillon

Sébastien Bonnefond

Camille Fabre

Sybille Roux
email: camille.roux@imag.fr

Caffiau

S Caffiau

C Roux

The SWEET-HOME Project: Audio Technology in Smart Homes to improve Well-being and Reliance

published or not. The documents may come

I. INTRODUCTION

Demographic change and ageing in developed countries imply challenges for the society to continue to improve the well being of its elderly and frail inhabitants. Since the dramatic evolution of Information and Communication Technologies (ICT), one way to achieve this aim is to promote the development of smart homes. In the health domain, a health smart home is a habitation equipped with a set of sensors, actuators, automated devices and centralised software controllers specifically designed for daily living task support, early detection of distress situations, remote monitoring and promotion of safety and well-being [START_REF] Chan | Smart homes -Current features and future perspectives[END_REF]. Among all the interaction and sensing technologies used in smart homes (e.g., infra-red sensors, contact doors, video cameras, RFID tags, etc.), audio processing technology has a great potential to become one of the major interaction modalities. Voice interfaces can be much more adapted to disabled people and the ageing population who have difficulties in moving or seeing than tactile interfaces (e.g., remote control) which require physical and visual interaction. Moreover, audio processing is particularly suited to distress situations [START_REF] Vacher | Development of audio sensing technology for ambient assisted living: Applications and challenges[END_REF].

The SWEET-HOME project which started in 2010 to address several of these issues is shortly introduced in Section II. The first task of the project was to acquire a large sample of situations for training and analysis purpose. Then, the core audio technologies were developed, namely, every day sound detection and distant speech recognition and incorporated into the real-time framework PATSH ("Plateforme d'Accueil de Traitement SWEET-HOME) which is presented in Section III together with the intelligent controller designed to interpret the voice command and make decision on the smart environment. The evaluation and the results are presented in Section IV before conclusion in Section V.

II. THE SWEET-HOME PROJECT

The SWEET-HOME project (sweet-home.imag.fr) is a French national supported research project aiming at designing a new smart home system based on audio technology focusing on three main aspects: to provide assistance via natural man-machine interaction (voice and tactile command), to ease social inclusion and to provide security reassurance by detecting situations of distress. If these aims are achieved, then the person will be able to pilot, from anywhere in the house, their environment at any time in the most natural way possible. The targeted users are elderly people who are frail but still autonomous. The architecture of the system is depicted in Figure 1. The input is composed of the information from the home automation network and information from 7 microphones (M1, M2,. . . , M7) transmitted through radio frequency channels. The audio signals are analysed by the PATSH framework described in Section III. Thus, information can be provided directly by the user (e.g., voice order) or via environmental sensors (e.g., temperature). The information coming from the home automation system is transmitted on-line to the intelligent controller (through several preprocessing steps). This controller captures all streams of data, interprets them and executes the required actions. the controller can send alert or information messages to a speech synthesiser in case of emergency or on request of the user ("Beware, the front door is not locked", "The temperature is 21 degrees").

Intelligent Controller

III. AUDIO ANALYSIS AND INTELLIGENT

CONTROLLER This section describes the audio analysis system and the intelligent controller.

A. Sound and speech analysis

For real-time audio analysis the PATSH framework (see Figure 1) was developed to manage the acquisition/processing flow of the sound events detected. In PATSH, several processors are plugged in and synchronised, such as the sound and speech analysis modules and the multichannel data acquisition card (NI-DAQ6220E). The system is composed of the following 4 stages:

1) the Sound Acquisition and Detection stage, detecting sound events (daily sounds or speech) from input streams; 2) the Sound/Speech Discrimination stage, discriminating speech from other sounds to extract voice commands; 3) the Sound Classification stage, recognizing daily living sounds; and 4) the Automatic Speech Recognition (ASR) stage, applying speech recognition to events classified as speech. The Sound/Speech Discrimination stage has a very important role: firstly, vocal orders must not be missed, secondly, daily living sounds must not be sent to the ASR because undesirable sentences could be recognized. In order to recognize only vocal orders and not all sentences uttered in the flat, all sentences shorter than 150 ms and longer than 2.2 seconds are ignored. These values were choosen after a statistical study on our data bases. In the same manner, the Signal to Noise Ratio (SNR) must be greater than 0 dB. For a complete description of the system the reader is referred to [START_REF] Sehili | Sound Environnment Analysis in Smart Home[END_REF].

B. Intelligent controller

The intelligent controller must not only analyse the recognized sentences but also monitor the smart home to detect risky situation. Its architecture and functioning are typical of the context-aware system model in pervasive environment.

The knowledge of the controller is defined using two semantic layers: the low-level ontology, devoted to the representation of raw data and network information description (e.g., state of switches and actuators); and the high-level ontology which represents concepts being used at the reasoning level. This separation between low and high levels makes it possible a high re-usability of the reasoning layer when the sensor network and the home have to be adapted [START_REF] Klein | Ontology-centred design of an ambient middleware for assisted living: The case of soprano[END_REF]. Outputs of the system include home automation orders but also interactions with the user in case of alert messages (e.g. warn the person about an unclosed front door).

The decision making module is based on an influence diagram extended to handle uncertain information and to take its decision based on the context. This decision is implemented through a Markov Logic Network which include the utility measure [START_REF] Nath | A language for relational decision theory[END_REF] (i.e., to which extend a decision is relevant given the situation) and a post-processing strategy to balance its decision given the uncertainty of the context. A detailed explanation about the decision method can be found in [START_REF] Chahuara | Context aware decision system in a smart home : knowledge representation and decision making using uncertain contextual information[END_REF]. The context provision is carried out through the collaboration of several processors executed under a service oriented approach, each one being specialized in a certain context aspect, such as location detection [START_REF]Location of an Inhabitant for Domotic Assistance Through Fusion of Audio and Non-Visual Data[END_REF] or activity recognition [START_REF] Chahuara | Using Markov Logic Network for On-line Activity Recognition from Non-Visual Home Automation Sensors[END_REF]. This data and knowledge centred approach ensures that all the processors are using the same data structure and that the meaning of each piece of information is clearly defined among all of them. In addition, the chosen architecture is more flexible than a classical pipeline of processors, making possible the easy insertion of new processors.

The decision module has two modes of action: 1) reactive when a vocal command is uttered and transmitted by PATSH, and 2) pro-active, when a risky situation is detected. The proactive mode is performed through situation recognition which triggers the decision module. The situation recognition is modelled by rules of the Event Condition Action (ECA) type implemented directly into the high-level ontology as Semantic Web Rule Language (SWRL) rules [START_REF] Horrocks | SWRL: A semantic web rule language combining OWL and RuleML[END_REF]. A typical example of such rule is showed below:

Situation. Event: The person open the front door to go out. Condition: one window is not closed. Action: Send a warning message using voice message.

SWRL rule DeviceEvent(?d), has associated object(?d, door), takes place in(?d, kitchen),state value(?d, open), Window(?w), located in(?w, bedroom),Application(?a), has application(?w,?a),curret state(?a,on) swlrb:moreThan(sqwrl:count(?w), 1) → current state(BedroomWindowsOpen, detected)

C. Vocal orders grammar

Our previous user study showed that targeted users prefer precise short sentences over more natural long sentences [START_REF] Portet | Design and evaluation of a smart home voice interface for the elderly: acceptability and objection aspects[END_REF]. As shown in Figure 2, each order belongs to one of three categories: initiate command, stop command and emergency call. Except for the emergency call, all command starts with a unique key-word that permits to know whether the person is talking to the smart home or not. In the following, we will use 'Nestor' as key-word: was estimated from all the sentences of the grammar. The recognized sentences are filterer at the input of the smart controller to determine their status, sentences not being in concordance with the grammar are rejected.

IV. EVALUATION The system was tested in a smart home during an experiment involving 13 persons. This section describes the pervasive environment, the experimental set-up, the experiment and the preliminary results.

A. Materials

Experiments were run in the DOMUS smart home that was designed by the Laboratory of Informatics of Grenoble [START_REF] Gallissot | Une plate-forme usage pour l'intégration de l'informatique ambiante dans l'habitat : Domus[END_REF]. Figure 3 shows the details of the flat. It is a thirty square meters suite flat including a bathroom, a kitchen, a bedroom and a study, all equipped with sensors and effectors. The flat is fully usable and can accommodate a dweller for several days. The technical architecture of DOMUS is based on the KNX bus system (standard ISO/IEC 14543). More than 150 sensors, actuators and information providers are managed in the flat (e.g., lighting, shutters, security systems, energy management, heating, etc.). A residential gateway architecture has been designed, supported by a virtual KNX layer seen as an Open Services Gateway Initiative service (OSGI) to guarantee the interoperability of the data coming and to allow the communication with virtual applications, such as activity tracking. The flat has also been equipped with 7 radio microphones for audio analysis with PATSH. A specialized communication device e-lio (www. technosens.fr/) from the Technosens company is used to initiate a communication with a relative.

B. Protocol

To validate the system in realistic conditions, we built scenarios in which every participant were asked to perform, at least once, some daily living activities chosen among the index of independence in Activities of Daily Living (ADL) [START_REF] Katz | Assessing self-maintenance: Activities of daily living, mobility, and instrumental activities of daily living[END_REF]. These activities include: (1) Sleeping; [START_REF] Vacher | Development of audio sensing technology for ambient assisted living: Applications and challenges[END_REF] Resting: listening to the radio; (3) Feeding: realizing and having a meal; and (4) Communicating: going out of the flat to do shopping and chatting with a relative thanks to e-lio. Therefore, this experiment allowed us to process realistic and representative audio events in conditions which are directly linked to usual daily living activities. Moreover, to evaluate the decision making, some specific situations were planned in the scenarios. For instance, for the decision regarding the activation of the light, given that the bedroom contains two lights (the ceiling and the bedside one) as well as the kitchen (above the dinning table and above the sink), the four following situations were planed:

1) Situation 1. The person is having a meal on the kitchen table. The most appropriate light is the one above the table. 2) Situation 2. The person is cleaning up the bedroom.

The most appropriate light is the ceiling one. 3) Situation 3. The person is cleaning the sink and doing the dishes. The most appropriate light is the one above the sink. 4) Situation 4. The person has just finished a nap. The most appropriate light is the bedside one. The participant had to use vocal orders to make the light on or off, open or close blinds, ask about temperature and ask to call his or her relative. The instruction was given to the participant to repeat the order up to 3 times in case of failure. If after 3 times the system did not react, the activation was performed remotely by Wizard of Oz.

C. Experiments

Thirteen healthy participants (including 4 women) were asked to perform the scenarios without condition on the duration. As in our previous experimentations [START_REF] Fleury | A french corpus of audio and multimodal interactions in a health smart home[END_REF], a visit, before the experiment, was organized to ensure that the participants will find all the items necessary to perform the seven ADLs. It was necessary to explain the right way to utter vocal orders and to use the e-lio system, too. The average age was 35 years (19-62, min-max) and the experiment lasted between 23mn 49s and 46mn 51s (total for all participants = 7h 1mn 42s). The scenario includes at least 15 vocal orders for each participant but more sentences were uttered because the person was asked to repeat the order in case of malfunctioning.

D. Results

Regarding the results of the global audio analysis system in the case of the participant S02, 73 events were classified as speech (including 6 sounds wrongly classified as speech), 10 speech events were classified as sounds. Among these sentences, 32 were vocal orders in accordance with the grammar (some were duplicated because recorded simultaneously by more than one microphone), 2 were bad recognized and only 2 were missed (not detected or classified as sound). Other sentences were system messages or part of the talk with the relative through the e-lio devices, a part of these sentences were too long and not analysed. Some sentences were uttered in presence of noise (curtain motor or radio).

Regarding the decision, it was highly dependent on the activity recognition. Indeed, the correct inference of the current activity is necessary to choose, e.g., the most appropriate light and its intensity. Among the 52 instances of the activity, 13 were sleeping, 13 were having a meal and 26 were cleaning. The accuracy of the activity recognition was 63.4% on this data set using models that consider 7 activities in total (meal, cleaning up, dressing, sleeping, leisure, hygiene, and communication) [START_REF] Chahuara | Using Markov Logic Network for On-line Activity Recognition from Non-Visual Home Automation Sensors[END_REF]. This is an interesting score given that the activities were not learned using the same data.

Table I show the accuracy of the decision making in each of the four situations. The first column 'with uncertainty' present the accuracy of the decision when taking the probability of all the activities into account to weight the utility of each decision in each situation in order to compute the expected utility. The second column 'no uncertainty' consisted in taking only the most probable activity into account. The lowest performance is in Situation 1. This can be explained by the fact that meal activity is often classified as clean-up. By contrast, the clean-up activity has a good recall mostly in the bedroom and thus leads to a correct decision. The results with and without uncertainty are similar as the cases in which the decision is different are only 4 out of 52. However, in one instance of Situation 3 the 'with uncertainty' case made it possible to find the correct solution. In that case, the highest activity probabilities were P (A = hygiene) = 0.2, P (A = dressing) = 0.16, P (A = sleeping) = 0.28, P (A = meal) = 0.06, P (A = communication) = 0.05 and P (A = leisure) = 0.17 whereas the ground truth was P (A = cleanup) = 0.08 in the kitchen. In this highly ambiguous situation, a correct decision was made because the most probable activities were the one with the highest degree of agitation. This helped the decision to choose the light with the highest intensity despite the most probable activity was sleeping.

V. CONCLUSIONS AND FUTURE WORKS

This paper presents an overview of the SWEET-HOME project which aims at designing a new smart home system based on audio technology. Several steps have been completed to provide real-time detection of sound of everyday living and speech recognition in the house. These modules were integrated inside the PATSH software to work in realtime. An intelligent controller makes the interface between the home automation network and PATSH in order to decide about the best action to make to adapt the smart home to the current situation. The experiment made in the smart home to evaluate the system showed promising results and validate the approach. This technology can benefit both the disabled and the elderly population that have difficulties in moving or seeing and want security reassurance.

Next steps in the project include the improvement of the current audio processing algorithms in order to detect in realtime the same vocal order recorded on several channels. This would improve the time response of the system in decreasing the number of sounds to analyse. Moreover, the time of analysis of the ASR is around 4s for each vocal order, which is far too long and biased the experiment. We are working in the optimisation of the ASR stage to perform new experiment with the targeted population with minimal technical bias.

Fig. 1 .

 1 Fig. 1. System architecture.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Excerpt of the grammar of the voice orders

TABLE I ACCURACY

 I OF THE DECISION MAKING.

	Accuracy:	with uncertainty	without uncertainty
	Situation 1	46%	54%
	Situation 2	100%	100%
	Situation 3	84%	69%
	Situation 4	54%	62%
	Total	71%	71%

VI. ACKNOWLEDGMENTS

The authors would like to thank the participants who took part to the different experiments. Thanks are also extended to D. Guérin to his support for PATSH development and to S. Pons for her support during the experiments inside the smart home.

This work is supported by the French National Research Agency (Agence Nationale de la Recherche /


D Nébouy

M Hébert

T Fournel

J.-L Lesur

Visual quality of printed surfaces: Study of homogeneity

Keywords: Print quality assessment, Halftoning, Quality index, Quality attributes, Uniformity, Gray-level co-occurrence matrix, Human Visual System, Appearance, Scanner, Printing technologies

This paper introduces a homogeneity assessment method for the printed versions of uniform color images. This parameter has been specifically selected as one of the relevant attributes of printing quality. The method relies on image processing algorithms from a scanned image of the printed surface, especially the computation of gray level cooccurrence matrices and of objective homogeneity attribute inspired of Haralick's parameters. The viewing distance is also taken into account when computing the homogeneity index. Resizing and filtering of the scanned image are performed in order to keep the level of details visible by a standard human observer at short and long distances. The combination of the obtained homogeneity scores on both high and low resolution images provides a homogeneity index, which can be computed for any printed version of a uniform digital image. We tested the method on several hardcopies of a same image, and compared the scores to the empirical evaluations carried out by non-expert observers who were asked to sort the samples and to place them on a metric scale. Our experiments show a good matching between the sorting by the observers and the score computed by our algorithm.

INTRODUCTION

Printing has evolved significantly the last few years, due to first a diversification of processes, and also to the digital revolution [START_REF] Kipphan | Handbook of Print Media[END_REF] including an automation of the processes without the expertise of the printman at the printing time. With this rapid evolution, the quality of color images reproduced on hard supports has considerably increased, and the amazing diversification of these techniques permitted by the new digital technologies provides a large panel of solutions for image reproduction. This leads to an issue for comparing the different techniques in terms of quality, an intuitive concept which is generally well appreciable visually but difficult to assess with objective, measurable values. Scientific concern for printing quality is recent: In [START_REF] Pedersen | Attributes of a New Image Quality Model for Color Prints[END_REF], the printing quality is explored by defining various quality attributes: color, lightness, sharpness, contrast, physical, artifacts. These attributes were estimated from psycho-visual experiments carried out by a panel of 15 observers with various expertise. In contrast with Pedersen et al. [START_REF] Pedersen | Framework for the Evaluation of Color Prints Using Image Quality Metrics[END_REF], [START_REF] Gong | Spatial pooling for measuring color printing quality attributes[END_REF], we believe that an important distinction should be made between printing quality (related to the degradation of the original digital image when transferred to the paper or plastic support) and image quality (including printing quality plus the perceived quality of the original image). Our opinion is that, contrast, lightness, and color rather qualify the image quality. The present work focuses on the homogeneity of tones i.e. the homogeneity of any surfaces on which is printed an originally uniform digital color image. This attribute is considered as relevant in the standard ISO 13660: 2001, the only international standard describing a wide range of attributes of image quality dedicated to quality assessment for binary or monochrome printing systems. [START_REF] Standard | IEC 13660 Information technology-Office equipment-Measurement of image quality attributes for[END_REF]. In this standard, mottling is defined as a type of homogeneity default of paper-like surfaces and an index is suggested in the spatial domain. A pass-band method was proposed [START_REF] Sadovnikov | Mottling Assessment of Solid Printed Areas and Its Correlation to Perceived Uniformity[END_REF] as an alternative able to better take into account the sensitivity of the human visual system. A pass-band assessment funded on wavelets, followed by a second-order statistics calculation on the gray-level co-occurrence matrix was achieved in [START_REF] Dube | Wavelet Analysis of Print Mottle[END_REF] to form an homogeneity index. Our aim here is to integrate human vision data into both the filtering step and the second-order statistics calculation. As in [START_REF] Sadovnikov | Mottling Assessment of Solid Printed Areas and Its Correlation to Perceived Uniformity[END_REF], our approach first consists in the digitalization of the printed surface using a high resolution flatbed office scanner, the scanned image being then processed in order to compute an homogeneity index. It then consists in comparing this score for many samples (several versions of a same digital image printed on different supports with different printing techniques) with the empirical evaluation by observers of these samples on a scale.

In the next sections, we first explain how the surfaces are digitalized and which tool of image processing we use before developing in details our method and its tunable parameters. We finally compare the score provided by the method with empirical assessment by observers on four sets of samples.

DIGITALIZATION OF THE PRINTS

For the digitalization of the samples, a scanner is used. An initial study on scanners have led to the conclusion that it is a reliable tool for image acquisition due to its remarquable reproductability, mainly performed by a uniform illumination, and good geometry that avoids specular reflections when the surface is flat. However, spatial and color distortions induced by the acquisition systems should be taken into account when assessing quality.

The color gamut of our human visual system is much larger than that of the scanner, which is unknown but can be assumed to be quite close to sRGB. Because of this, the color distortion problem is twofold: First, the so-called "seen by the scanner" colors are displaced within the very gamut of the scanner, but also, it may appears that different printable colors are out of the gamut of the scanner and thus perceived by the scanner as having the same color. A color calibration would seem necessary to solve at least the first issue. However, we noticed that very often after color correction, images are noisier than before correction. As part of our work (the assessment of uniformity), we chose not to apply color correction to avoid the risk of disrupting our measure.

Regarding the spatial distortions, there is a significant horizontal drift. In the vertical direction, we do not notice geometric drift, the motor is properly position controlled. This drift can be corrected by adjusting the position of the pixels by polynomial interpolation or spline on each line of the captured image. But such treatment would create a change in the value of certain pixels on each line, the most affected by this drift being on the edges of each line, and would lead to a loss of the native resolution of the image as well as local color unwanted distortions. It is preferable to analyze the scanner profile derives and locate the sample on the scanner glass in an area where the drifts are minimal.

We recommend scanning samples at a high resolution (e.g. 2400 dpi), and then reducing to 1200 dpi, the initial resolution of our method, in order to perform controlled down-sampling.

BACKGROUND: IMAGE HOMOGENEITY ANALYSIS

The printed copy of a uniform color patch can be considered as a texture. In image processing, a classical tool in texture characterization relies on Haralick's parameters derived from the gray-level co-occurrence matrix (GLCM) also called gray-tone spatial-dependence matrix, first presented in [START_REF] Haralick | Textural Features for Image Classification[END_REF]. In contrast with first-order statistic tools (histogram, mean, standard deviation, skewness, and kurtosis), the GLCM is a statistic tool of 2 nd order [START_REF] Srinivasan | Statistical texture analysis[END_REF]. First order statistics are representative of gray-level distribution of the pixels in the image regardless of their spatial arrangement, whereas second-order statistics involve two pixels simultaneously.

The construction of the GLCM t H is defined as:

  0..2 , 0..2 (, ) nn bb ij hij         tt H (1)
with

    2 ( , ) , ; ( ) and ( ) hij C a r d pp I f p i f p j     t tt ( 2
)
where I is a rectangular image encoded on b n bits, f is an application from 0..2 1 b n II      , Card means the cardinality of a set, t is a given translation vector, and   , pp  t is a pair of pixels in image I . The vector t may be specified by a unit vector w and a distance D :

. D  tw [START_REF] Pedersen | Framework for the Evaluation of Color Prints Using Image Quality Metrics[END_REF] Starting from an initialization to the null matrix, t H is incrementally built by considering in turn every pair of pixels

 

, pp  t in image I , and incrementing by one the entry   () , ( ) hfpfp  t t (see Figure 1).

After normalization of this matrix, each entry   , hij t means the probability to transit, along a vector t , from a pixel with value i to a pixel with value j :

       , , hij pij Card I T I    t t t (4)
with T x means the translation operator of a vector x . On such GLCM, Haralick defines several texture attributes among which an homogeneity attribute comprised between 0 and 1, defined as:

  22 2 00 1 , 1( ) nn bb ij Hom p i j ij     t (5)
If the GLCM t H of image I is computed on a uniform image (top left image on Figure 2) or on a periodical image with period t (top middle image), then for every pair of pixels

      , fp fp  t , we have    fp fp  t . Therefore,
only entries in the diagonal of t H are incremented; and 1 Hom  . In the other cases (for example for the top right image of Figure 2),t H is not a diagonal matrix, and 01 Hom   .

Formula ( 5) may be replaced with other ones, defined so as to give different weight to entries more distant from the diagonal. For example, we may generalize Eq. ( 5) by defining the following function weighting the transition probability

  , t pij: , 1 :
1

ab b wx ax    (6)
where a and b are tunable real values (see Figure 3), and x ij   .

Using , ab w with different a and b values than the ones used in Eq. ( 5), we can tune the response of the homogeneity score to variations of the pixel values in the image.

OVERVIEW OF THE METHOD

The Homogeneity assessment method that we propose relies on the following steps:

 Preprocessing, including image resizing and filtering in order to simulate near vision at High Resolution (HR), and far vision at Low Resolution (LR),  Application of the Image Homogeneity Analysis method described in Section 3.

 Score computation from the GLCM on both HR and LR images.

PREPROCESSING

This section intends to modify the raw scanned image into a new one, closer to human perception. A three step preprocessing of the scanned images is performed: Creation of custom color channels from the initial RGB channels; filtering in order to take account of human visual resolution; quantization (optional).

For the HR image, a 1200 pixel per inch (ppi) image is recommended. A bicubic interpolation (scale of 1/2 if the input image is 2400 ppi) with antialiasing method is selected to perform image resizing.

The first step is to change the color space. Scanner yields R,G,B values for each pixel. We could process these R,G, and B channels each ones as a grayscale image, but they do not contain the optimal information. In order to be more

The second step corresponds to the simulation of the blur produced by our human eyes. According to Hermann von Helmholtz's observations [START_REF] Helmholtz | Helmholtz's Treatise on Physiological Optics[END_REF], we consider the definition of the "normal visual acuity" as the ability to resolve a spatial pattern separated by a visual angle α of one minute of arc. Thus, giving an observation distance d (in cm) and a scan resolution of the image s r (in ppi), we can define the convolution kernel of the blurring filter (Figure 4). Lastly, a quantization step (e.g. on 8

b n  bits for each channel), can eventually be performed, with the consequence to decrease the size of the GLCM but also to lose the least significant bits (i.e. small grey-level variations in the image).

For the low resolution images, we perform the same 3 steps with a 240 ppi resolution image (scale of 1/10 if the input is 2400 ppi), with 1.03   pixels and 7 t  , which corresponds to the image viewed by an observer at about 40 cm. The preprocessing step leads to 4 images, denoted as    ,,,

HRL RH RL R ab ab

YYhh .

APPLICATION OF THE IMAGE HOMOGENEITY ANALYSIS METHOD

We apply the method presented in Section 3 on the four images    ,,,

HRL RH RL R ab ab

YYhh obtained in the preprocessing step.

For each of the two HR images, we choose the following translation vectors 1 t and 2 t to compute respectively GLCMs We then compute a GLCM M as the average 1 M and 2 M :

2

     12 MM M ( 14
)
We also compute a chroma value on the HR image, denoted as   HR c and defined as the average of the chroma values computed from Eq. ( 7) on each pixel of the image. From matrix M , we compute the homogeneity attribute Hom1 according to Eq. ( 6) with a value of a denoted   For each of the two LR images, we perform similar processing as above with   35 LR D  pixels (corresponding to a physical displacement of about 0.37 cm on 240 ppi images), and compute the Hom2 attribute using the   ,4

CONSTRUCTION OF THE OBJECTIVE SCORE

The computation of homogeneity attribute on the GLCM of    ,,,

HRL RH RL R ab ab

YYhh yields four values to be combined in order to obtain the final homogeneity score. Let's denote these four values as

        ,,, HRH R L RL R ab ab YY hh HHHH .
Our experience shows that relevant homogeneity score in respect to human perception is given by the following empirical formula: computed in Section 6 because we noticed that the influence of the homogeneity parameter on the hue image is lower when the colors are not saturated. The weight attributed to the luminance channel is higher than the one attributed to the hue channel ( 01 c   , and in practice, on prints, it is rather 00 . 6 c  ) as it is commonly done in video compression for digital television video encoding where the luminance channel has a double bandwidth than the chromatic ones (See for example [START_REF] Itu | Studio encoding parameters of digital television for standard 4:3 and wide-screen 16:9 aspect ratios[END_REF]).

                min , 11

REFERENCED PARAMETERS OF THE METHOD

The referenced parameters of our method are:

 Quantization (images encoded on b n number of bits): The less is quantization, the more careful details are kept. It seems better to keep all the details. Instead of decreasing b n below the usual value of 8 bits, we prefer adjusting the a and b parameters in formula (6).

 Size of blurring kernel: it is set by the considered viewing distance d , which is about 15 cm for near vision, and 40 cm for far vision.

 Scan resolution: The size of one pixel in HR images is about 21 μm after downsampling original 2400 ppi scanned image to 1200 ppi. This appears sufficiently accurate in comparison to the human vision and to the finest commercial printing systems available today. t : horizontal and vertical vectors used to compute the GLCM. These two directions have been selected because they coincide to the orientation of classical defects of most printing systems, because defaults often occur in the printing direction (vertical) and in the perpendicular direction of it.

EXPERIMENTAL TESTING AND COMPARISON WITH PSYCHO-VISUAL ASSESSMENT

In order to verify the relevance of the proposed score in respect to human observer evaluation, we carried out the following experiment: Uniform color patches have been printed with various printers (laser, inkjet, retransfer) and on various supports (office paper, glossy photograph paper, APCO paper, white polymer). For the digitalization step, all samples are scanned at 2400 ppi in TIFF format, 24 bits, by the EPSON Scan Ver. 3.81 FR, 2002 driver, with a EPSON Perfection V700 PHOTO scanner. The "No color correction" option is selected. An overview of some samples is presented below. Figure 6. An overview of some scanned samples (in lower resolution than the original ones). s_02 is printed on glossy photograph paper, s_01 and s_08 on office paper, s_10 and s_11 on APCO paper.

The assessors were submitted to a ranking test, using an ordinale scale, unstructured, as it is presented by François Sauvageot in [START_REF] Depledt | Evaluation sensorielle Manuel méthodologique[END_REF]. The assessors were presented randomly the samples to sort in respect to their perceived homogeneity, and were asked to place them on a one meter scale. The positions of the samples (in meter) are indicated by the dots in Figure 7, where one dot shape is attributed to each observer. Note that no rescaling of the positions has been performed: the dispersion of the positions for one sample partly comes from the different scaling adopted by the observers in their evaluation. The continuous line in Figure 7 represents the score computed with our method, represented by formula (15). We see that the scores follows fairly well the positioning of the samples by the observers on the metric scale. In order to assess the performance of our method, six different classes are defined by selecting six threshold values of the score (0.8, 0.7, 0.6, 0.52, 0.3, 0). Assessors scorings are also quantified into this six grades classification. The matching of the class given by observers and by our method is presented in Figure 8.

Figure 1 .Figure 2 .

 12 Figure 1. Scanning of the image pixels for every point   , p px y in order to build the gray level co-occurrence matrix.

2 w

 2 Figure 3. , ab w function for three different sets of   , ab values. Original Haralick's attribute (denoted as Hom) corresponds

Figure 4 .

 4 Figure 4. Definition of  knowing the visual angular resolution

1 uy ,  2 ux

 12 to a physical displacement of about 0.5 mm on 1200 ppi image), and  the horizontal, respectively vertical unit vectors (see Figure1).

  values of the LR image. These values for D , a , and b , have been chosen because they seem the most relevant in our study according to physical and empirical considerations.

  Figure 5. Variation of the   HR a

  chroma is taken into account in this formula through the parameters   HR c and   LR c

  must be representative of the characteristic size of the inhomogeneities that we wish to highlight, in respect to near sight and far sight observations.

Figure 7 .

 7 Figure 7. Comparison between objective and subjective evaluation for (a): a blue set of samples, (b): an orange set of samples, (c): a red set of samples, and (d): a gray set of samples.

Figure 8 .

 8 Figure 8. Percentage of samples classified in the same of six categories by assessors (subjective evaluation) and by our scoring method.

  consistent with human perception, we prefer converting the RGB values into CIELAB values by considering that the input RGB values are represented in the Adobe RGB (1998) color space, then converting into CIE1931XYZ for a 2° standard observer, and finally converting in CIELAB by selecting the D65 illuminant. We analyze the *

						L (lightness)
	channel, and the ab h (hue) channel. CIELAB color space is chosen because it is a color appearance model [10] (pages
	160 -161). The * a and * b Cartesian coordinates are transformed into cylindrical coordinates * ab C and ab h , respectively
	defined as
	** 2 ab Ca 	b	* 2	(7)
	ab h		tan	* * a b   1  	(8)
	Implementation remark: function is available in most computation softwares, often denoted as   atan y x function yields an angle in 22 ,        atan 2 , which is not defined if   y x , yielding an angle in  x  . Alternative 0  ,    , defined
	for every y and x .
	Note that it is also possible to define the luminance image directly from RGB values, using the NTSC luminance (Y)
	formula, first presented in Ref. [11]:
	0.299		0.587		0.114

YRGB   

CONCLUSION

The index computed according to the method that we propose for the assessment of the homogeneity of printed colors is in accordance with the empirical assessment by a small panel of non-expert observers. In this sense, if the trend is confirmed for a larger set of observers, the algorithm that we have developed will appear consistent with the human visual perception. This homogeneity index is a first attribute which should be combined with other attributes in order to get a global quality score for printed natural images. It may also help to compare the performance of different printing systems in terms of visual quality.