Bruno Scherrer
email: bruno.scherrer@inria.fr

Mohammad Ghavamzadeh
email: mohammad.ghavamzadeh@inria.fr

Victor Gabillon
email: victor.gabillon@inria.fr

Matthieu Geist
email: matthieu.geist@supelec.fr

Approximate Modified Policy Iteration

Modified policy iteration (MPI) is a dynamic programming (DP) algorithm that contains the two celebrated policy and value iteration methods. Despite its generality, MPI has not been thoroughly studied, especially its approximation form which is used when the state and/or action spaces are large or infinite. In this paper, we propose three implementations of approximate MPI (AMPI) that are extensions of well-known approximate DP algorithms: fitted-value iteration, fitted-Q iteration, and classification-based policy iteration. We provide error propagation analysis that unifies those for approximate policy and value iteration. For the classificationbased implementation, we develop a finitesample analysis that shows that MPI's main parameter allows to control the balance between the estimation error of the classifier and the overall value function approximation.

Introduction

Modified Policy Iteration (MPI) [START_REF] Puterman | Modified policy iteration algorithms for discounted Markov decision problems[END_REF] is an iterative algorithm to compute the optimal policy and value function of a Markov Decision Process (MDP). Starting from an arbitrary value function v 0 , it generates a sequence of value-policy pairs

π k+1 = G v k (greedy step) (1) v k+1 = (T π k+1 ) m v k (evaluation step) (2)
where G v k is a greedy policy w.r.t. v k , T π k is the Bellman operator associated to the policy π k , and m ≥ 1 is a parameter. MPI generalizes the well-known dynamic programming algorithms Value Iteration (VI) and Policy Iteration (PI) for values m = 1 and m = ∞, respectively. MPI has less computation per iteration than PI (in a way similar to VI), while enjoys the faster convergence of the PI algorithm [START_REF] Puterman | Modified policy iteration algorithms for discounted Markov decision problems[END_REF].

In problems with large state and/or action spaces, approximate versions of VI (AVI) and PI (API) have been the focus of a rich literature (see e.g., [START_REF] Bertsekas | Neuro-Dynamic Programming[END_REF][START_REF] Szepesvári | Reinforcement learning algorithms for MDPs[END_REF]. The aim of this paper is to show that, similarly to its exact form, approximate MPI (AMPI) may represent an interesting alternative to AVI and API algorithms.

In this paper, we propose three implementations of AMPI (Sec. 3) that generalize the AVI implementations of [START_REF] Ernst | Tree-based batch mode reinforcement learning[END_REF]; [START_REF] Antos | Fitted Qiteration in continuous action-space MDPs[END_REF]; [START_REF] Munos | Finite-time bounds for fitted value iteration[END_REF] and the classification-based API algorithm of [START_REF] Lagoudakis | Reinforcement learning as classification: Leveraging modern classifiers[END_REF]; [START_REF] Fern | Approximate policy iteration with a policy language bias: Solving relational markov decision processes[END_REF]; [START_REF] Lazaric | Analysis of a classification-based policy iteration algorithm[END_REF]; [START_REF] Gabillon | Classification-based policy iteration with a critic[END_REF]. We then provide an error propagation analysis of AMPI (Sec. 4), which shows how the L p -norm of its performance loss can be controlled by the error at each iteration of the algorithm. We show that the error propagation analysis of AMPI is more involved than that of AVI and API. This is due to the fact that neither the contraction nor monotonicity arguments, that the error propagation analysis of these two algorithms rely on, hold for AMPI. The analysis of this section unifies those for AVI and API and is applied to the AMPI implementations presented in Sec. 3. We detail the analysis of the classification-based implementation of MPI (CBMPI) of Sec. 3 by providing its finite sample analysis in Sec. 5. Our analysis indicates that the parameter m allows us to balance the estimation error of the classifier with the overall quality of the value approximation. We report some preliminary results of applying CBMPI to standard benchmark problems and comparing it with some existing algorithms in Appendix. G.

Background

We consider a discounted MDP �S, A, P, r, γ�, where S is a state space, A is a finite action space, P (ds � |s, a), for all (s, a), is a probability kernel on S, the reward function r : S × A → R is bounded by R max , and γ ∈ (0, 1) is a discount factor. A deterministic policy is defined as a mapping π : S → A. For a policy π, we may write

r π (s) = r � s, π(s) � and P π (ds � |s) = P � ds � |s, π(s) � .
The value of policy π in a state s is defined as the expected discounted sum of rewards received starting from state s and following the policy π, i.e., v π (s

) = E � � ∞ t=0 γ t r π (s t )| s 0 = s, s t+1 ∼ P π (•|s t ) � .
Similarly, the action-value function of a policy π at a state-action pair (s, a), Q π (s, a), is the expected discounted sum of rewards received starting from state s, taking action a, and then following the policy. Since the rewards are bounded by R max , the values and action-values should be bounded by

V max = Q max = R max /(1 -γ).
The Bellman operator T π of policy π takes a function f on S as input and returns the function

T π f defined as ∀s, [T π f ](s) = E � r π (s) + γf (s � ) | s � ∼ P π (.|s) �
, or in compact form, T π f = r π + γP π f . It is known that v π is the unique fixed-point of T π . Given a function f on S, we say that a policy π is greedy w.r.t. f , and write it as π = G f , if ∀s, (T π f )(s) = max a (T a f )(s), or equivalently T π f = max π � (T π � f ). We denote by v * the optimal value function. It is also known that v * is the unique fixed-point of the Bellman optimality operator

T : v → max π T π v = T G(v) v,
and that a policy π * that is greedy w.r.t. v * is optimal and its value satisfies v π * = v * .

Approximate MPI Algorithms

In this section, we describe three approximate MPI (AMPI) algorithms. These algorithms rely on a function space F to approximate value functions, and in the third algorithm, also on a policy space Π to represent greedy policies. In what follows, we describe the iteration k of these iterative algorithms.

AMPI-V

For the first and simplest AMPI algorithm presented in the paper, we assume that the values v k are represented in a function space F ⊆ R |S| . In any state s, the action π k+1 (s) that is greedy w.r.t. v k can be estimated as follows:

π k+1 (s) = arg max a∈A 1 M � M � j=1 r (j) a + γv k (s (j) a ) � , (3)
where ∀a ∈ A and 1 ≤ j ≤ M , r (j) a and s

(j) a are samples of rewards and next states when action a is taken in state s. Thus, approximating the greedy action in a state s requires M |A| samples. The algorithm works as follows. It first samples N states from a distribution µ, i.e., {s (i) } N i=1 ∼ µ. From each sampled state s (i) , it generates a rollout of size m, i.e., � s (i) , a

(i) 0 , r (i) 0 , s (i) 1 , . . . , a (i) m-1 , r (i) m-1 , s (i) m �
, where a (i) t is the action suggested by π k+1 in state s (i) t , computed using Eq. 3, and r (i) t and s (i) t+1 are the reward and next state induced by this choice of action. For each s (i) , we then compute a rollout estimate

� v k+1 (s (i) ) = � m-1 t=0 γ t r (i) t + γ m v k (s (i) m ), which is an un- biased estimate of �� T π k+1 � m v k � (s (i)
). Finally, v k+1 is computed as the best fit in F to these estimates, i.e.,

v k+1 = Fit F � � � s (i) , � v k+1 (s (i) ) � � N i=1 � .
Each iteration of AMPI-V requires N rollouts of size m, and in each rollout any of the |A| actions needs M samples to compute Eq. 3. This gives a total of N m(M |A|+1) transition samples. Note that the fitted value iteration algorithm [START_REF] Munos | Finite-time bounds for fitted value iteration[END_REF]) is a special case of AMPI-V when m = 1.

AMPI-Q

In AMPI-Q, we replace the value function v : S → R with an action-value function Q : S × A → R. The Bellman operator for a policy π at a state-action pair (s, a) can then be written as

[T π Q](s, a) = E � r π (s, a) + γQ(s � , π(s � ))|s � ∼ P (•|s, a) � ,
and the greedy operator is defined as

π = G Q ⇔ ∀s π(s) = arg max a∈A Q(s, a).
In AMPI-Q, action-value functions Q k are represented in a function space F ⊆ R |S×A| , and the greedy action w.r.t. Q k at a state s, i.e., π k+1 (s), is computed as

π k+1 (s) ∈ arg max a∈A Q k (s, a). (4)
The evaluation step is similar to that of AMPI-V, with the difference that now we work with stateaction pairs. We sample N state-action pairs from a distribution µ on S × A and build a rollout set

Construct the rollout set D k = {s (i) } n i=1 , s (i) iid ∼ µ for all states s (i) ∈ D k do Perform a rollout and return � v k (s (i) ) end for Construct the rollout set D � k = {s (i) } N i=1 , s (i) iid ∼ µ for all states s (i) ∈ D �
k and actions a ∈ A do for j = 1 to M do Perform a rollout and return

R j k (s (i) , a) end for � Q k (s (i) , a) = 1 M � M j=1 R j k (s (i) , a) end for • Approximate value function: v k ∈ argmin v∈F � L F k (� µ; v) (regression)
• Approximate greedy policy:

π k+1 ∈ argmin π∈Π � L Π k (� µ; π) (classification)
end for

Figure 1. The pseudo-code of the CBMPI algorithm.

D k = {(s (i) , a (i) )} N i=1 , (s (i) , a (i) ) ∼ µ. For each (s (i) , a (i) ) ∈ D k , we generate a rollout of size m, i.e., � s (i) , a (i) , r (i) 0 , s (i)
1 , a

(i) 1 , • • • , s (i) m , a (i) m
� , where the first action is a (i) , a (i) t for t ≥ 1 is the action suggested by π k+1 in state s (i) t computed using Eq. 4, and r (i) t and s (i) t+1 are the reward and next state induced by this choice of action. For each (s (i) , a (i) ) ∈ D k , we then compute the rollout estimate

� Q k+1 (s (i) , a (i) ) = m-1 � t=0 γ t r (i) t + γ m Q k (s (i) m , a (i) m ), which is an unbiased estimate of � (T π k+1 ) m Q k � (s (i) , a (i) ). Finally, Q k+1 is the best fit to these estimates in F, i.e., Q k+1 = Fit F � � � (s (i) , a (i) ), � Q k+1 (s (i) , a (i) ) � � N i=1 � .
Each iteration of AMPI-Q requires N m samples, which is less than that for AMPI-V. However, it uses a hypothesis space on state-action pairs instead of states. Note that the fitted-Q iteration algorithm [START_REF] Ernst | Tree-based batch mode reinforcement learning[END_REF][START_REF] Antos | Fitted Qiteration in continuous action-space MDPs[END_REF] is a special case of AMPI-Q when m = 1.

Classification-Based MPI

The third AMPI algorithm presented in this paper, called classification-based MPI (CBMPI), uses an ex-plicit representation for the policies π k , in addition to the one used for value functions v k . The idea is similar to the classification-based PI algorithms [START_REF] Lagoudakis | Reinforcement learning as classification: Leveraging modern classifiers[END_REF][START_REF] Fern | Approximate policy iteration with a policy language bias: Solving relational markov decision processes[END_REF][START_REF] Lazaric | Analysis of a classification-based policy iteration algorithm[END_REF][START_REF] Gabillon | Classification-based policy iteration with a critic[END_REF] in which we search for the greedy policy in a policy space Π (defined by a classifier) instead of computing it from the estimated value or action-value function (like in AMPI-V and AMPI-Q).

In order to describe CBMPI, we first rewrite the MPI formulation (Eqs. 1 and 2) as

v k = (T π k ) m v k-1 (evaluation step) (5) π k+1 = G � (T π k ) m v k-1 � (greedy step) (6)
Note that in the new formulation both v k and π k+1 are functions of (T π k ) m v k-1 . CBMPI is an approximate version of this new formulation. As described in Fig. 1, CBMPI begins with arbitrary initial policy

π 1 ∈ Π and value function v 0 ∈ F. 1 At each iteration k, a new value function v k is built as the best approx- imation of the m-step Bellman operator (T π k ) m v k-1 in F (evaluation step). This is done by solving a re- gression problem whose target function is (T π k ) m v k-1 .
To set up the regression problem, we build a rollout set D k by sampling n states i.i.d. from a distribution µ.2 For each state

s (i) ∈ D k , we generate a roll- out � s (i) , a (i) 0 , r (i) 0 , s (i) 1 , . . . , a (i)
m-1 , r

(i) m-1 , s (i) m � of size m, where a (i) t = π k (s (i) t ),
v k (s (i) ) of � (T π k ) m v k-1 � (s (i) ) as � v k (s (i) ) = m-1 � t=0 γ t r (i) t + γ m v k-1 (s (i) m ), (7)
and use it to build a training set ��

s (i) , � v k (s (i) ) �� n i=1
. This training set is then used by the regressor to compute v k as an estimate of (T π k ) m v k-1 .

The greedy step at iteration k computes the policy π k+1 as the best approximation of G

� (T π k ) m v k-1 � by solv- ing a cost-sensitive classification problem. From the definition of a greedy policy, if π = G � (T π k ) m v k-1 � , for each s ∈ S, we have � T π (T π k ) m v k-1 � (s) = max a∈A � T a (T π k ) m v k-1 � (s). (8) By defining Q k (s, a) = � T a (T π k ) m v k-1 � (s)
, we may rewrite Eq. 8 as

Q k � s, π(s) � = max a∈A Q k (s, a). (9)
The cost-sensitive error function used by CBMPI is of the form

L Π π k ,v k-1 (µ; π) = � S � max a∈A Q k (s, a) -Q k � s, π(s) � � µ(ds).
To simplify the notation we use

L Π k instead of L Π π k ,v k-1 .
To set up this cost-sensitive classification problem, we build a rollout set D � k by sampling N states i.i.d. from a distribution µ. For each state s (i) ∈ D � k and each action a ∈ A, we build M independent rollouts of size m + 1, i.e.,3 � s (i) , a, r

(i,j) 0 , s (i,j) 1 , a (i,j) 1 , . . . , a (i,j) m , r (i,j) m , s (i,j) m+1 � M j=1 ,
where for t ≥ 1, a

(i,j) t = π k (s (i,j) t
), and r (i,j) t and s (i,j) t+1 are the reward and next state induced by this choice of action. From these rollouts, we compute an unbiased estimate of

Q k (s (i) , a) as � Q k (s (i) , a) = 1 M � M j=1 R j k (s (i) , a)
where

R j k (s (i) , a) = m � t=0 γ t r (i,j) t + γ m+1 v k-1 (s (i,j) m+1 ).
Given the outcome of the rollouts, CBMPI uses a costsensitive classifier to return a policy π k+1 that minimizes the following empirical error

� L Π k (� µ; π) = 1 N N � i=1 � max a∈A � Q k (s (i) , a) -� Q k � s (i) , π(s (i) ) � � ,
with the goal of minimizing the true error L Π k (µ; π). Each iteration of CBMPI requires nm+M |A|N (m+1) (or M |A|N (m + 1) in case we reuse the rollouts, see Footnote 3) transition samples. Note that when m tends to ∞, we recover the DPI algorithm proposed and analyzed by [START_REF] Lazaric | Analysis of a classification-based policy iteration algorithm[END_REF].

Error propagation

In this section, we derive a general formulation for propagation of error through the iterations of an AMPI algorithm. The line of analysis for error propagation is different in VI and PI algorithms. VI analysis is based on the fact that this algorithm computes the fixed point of the Bellman optimality operator, and this operator is a γ-contraction in max-norm [START_REF] Bertsekas | Neuro-Dynamic Programming[END_REF][START_REF] Munos | Performance bounds in Lp-norm for approximate value iteration[END_REF]. On the other hand, it can be shown that the operator by which PI updates the value from one iteration to the next is not a contraction in max-norm in general. Unfortunately, we can show that the same property holds for MPI when it does not reduce to VI (i.e., m > 1).

Proposition 1. If m > 1, there exists no norm for which the operator that MPI uses to update the values from one iteration to the next is a contraction.

Proof. Consider a deterministic MDP with two states {s 1 , s 2 }, two actions {change, stay}, rewards r(s 1 ) = 0, r(s 2 ) = 1, and transitions P ch (s 2 |s 1 ) = P ch (s 1 |s 2 ) = P st (s 1 |s 1 ) = P st (s 2 |s 2 ) = 1. Consider the following two value functions v = (�, 0) and v � = (0, �) with � > 0. Their corresponding greedy policies are π = (st, ch) and π � = (ch, st), and the next iterates of v and v � can be computed as

(T π ) m v = � γ m � 1 + γ m � � and (T π � ) m v � = � γ-γ m 1-γ + γ m � 1-γ m 1-γ + γ m � � . Thus, (T π � ) m v � -(T π ) m v = � γ-γ m 1-γ γ-γ m 1-γ � while v � -v = � -� � � . Since � can be arbitrarily small, the norm of (T π � ) m v � -(T π ) m v can be arbitrarily larger than the norm of v -v � as long as m > 1.
We also know that the analysis of PI usually relies on the fact that the sequence of the generated values is non-decreasing [START_REF] Bertsekas | Neuro-Dynamic Programming[END_REF][START_REF] Munos | Error bounds for approximate policy iteration[END_REF]. Unfortunately, it can be easily shown that for m finite, the value functions generated by MPI may decrease (it suffices to take a very high initial value).

It can be seen from what we just described and Proposition 1 that for m � = 1 and ∞, MPI is neither contracting nor non-decreasing, and thus, a new line of proof is needed for the propagation of error in this algorithm.

To study error propagation in AMPI, we introduce an abstract algorithmic model that accounts for potential errors. AMPI starts with an arbitrary value v 0 and at each iteration k ≥ 1 computes the greedy policy w.r.t. v k-1 with some error � � k , called the greedy step error. Thus, we write the new policy π k as

π k = � G � � k v k-1 . (10)
Eq. 10 means that for any policy π � ,

T π � v k-1 ≤ T π k v k-1 + � � k .
AMPI then generates the new value function v k with some error � k , called the evaluation step error

v k = (T π k ) m v k-1 + � k . ( 11
)
Before showing how these two errors are propagated through the iterations of AMPI, let us first define them in the context of each of the algorithms presented in Section 3 separately.

AMPI-V: � k is the error in fitting the value function v k . This error can be further decomposed into two parts: the one related to the approximation power of F and the one due to the finite number of samples/rollouts. � � k is the error due to using a finite number of samples M for estimating the greedy actions.

AMPI-Q: � � k = 0 and � k is the error in fitting the state-action value function Q k .

CBMPI: This algorithm iterates as follows:

v k = (T π k ) m v k-1 + � k π k+1 = � G � � k+1 [(T π k ) m v k-1
] Unfortunately, this does not exactly match with the model described in Eqs. 10 and 11. By introducing the auxiliary variable

w k Δ = (T π k ) m v k-1 , we have v k = w k + � k ,
and thus, we may write

π k+1 = � G � � k+1 [w k ] . ( 12
) Using v k-1 = w k-1 + � k-1 , we have w k = (T π k ) m v k-1 = (T π k ) m (w k-1 + � k-1 ) = (T π k ) m w k-1 + (γP π k ) m � k-1 . (13)
Now, Eqs. 12 and 13 exactly match Eqs. 10 and 11 by replacing v k with w k and � k with (γP

π k ) m � k-1 .
The rest of this section is devoted to show how the errors � k and � � k propagate through the iterations of an AMPI algorithm. We only outline the main arguments that will lead to the performance bound of Thm. 1 and report most proofs in [START_REF] Scherrer | Approximate modified policy iteration[END_REF]. We follow the line of analysis developped by [START_REF] Thiery | Performance bound for approximate optimistic policy iteration[END_REF]. The results are obtained using the following three quantities:

1) The distance between the optimal value function and the value before approximation at the k th iteration:

d k Δ = v * -(T π k ) m v k-1 = v * -(v k -� k ).
2) The shift between the value before approximation and the value of the policy at the k th iteration:

s k Δ = (T π k ) m v k-1 -v π k = (v k -� k ) -v π k .
3) The Bellman residual at the k th iteration:

b k Δ = v k -T π k+1 v k .
We are interested in finding an upper bound on the

loss l k Δ = v * -v π k = d k + s k .
To do so, we will upper bound d k and s k , which requires a bound on the Bellman residual b k . More precisely, the core of our analysis is to prove the following point-wise inequalities for our three quantities of interest.

Lemma 1 (Proof in Appendix. A). Let k ≥ 1, x k Δ = (I -γP π k )� k + � � k+1 and y k Δ = -γP π * � k + � � k+1 . We have: b k ≤ (γP π k ) m b k-1 + x k , d k+1 ≤ γP π * d k + y k + m-1 � j=1 (γP π k+1 ) j b k , s k = (γP π k ) m (I -γP π k ) -1 b k-1 .
Since the stochastic kernels are non-negative, the bounds in Lemma 1 indicate that the loss l k will be bounded if the errors � k and � � k are controlled. In fact, if we define � as a uniform upper-bound on the errors

|� k | and |� � k |, the first inequality in Lemma 1 implies that b k ≤ O(�)
, and as a result, the second and third inequalities give us d k ≤ O(�) and s k ≤ O(�). This means that the loss will also satisfy l k ≤ O(�).

Our bound for the loss l k is the result of careful expansion and combination of the three inequalities in Lemma 1. Before we state this result, we introduce some notations that will ease our formulation.

Definition 1. For a positive integer n, we define P n as the set of transition kernels that are defined as follows:

1) for any set of n policies {π 1 , . . . , π n }, (γP π1 )(γP π2 ) . . . (γP πn ) ∈ P n , 2) for any α ∈ (0, 1) and (P 1 , P 2 ) ∈ P n × P n , αP 1 + (1α)P 2 ∈ P n . Furthermore, we use the somewhat abusive notation Γ n for denoting any element of P n . For example, if we write a transition kernel P as

P = α 1 Γ i + α 2 Γ j Γ k = α 1 Γ i +α 2 Γ j+k ,
it should be read as there exist P 1 ∈ P i , P 2 ∈ P j , P 3 ∈ P k , and P 4 ∈ P k+j such that P = α 1 P 1 + α 2 P 2 P 3 = α 1 P 1 + α 2 P 4 .

Using the notation introduced in Definition 1, we now derive a point-wise bound on the loss.

Lemma 2 (Proof in Appendix. B). After k iterations, the losses of AMPI-V and AMPI-Q satisfy

l k ≤ 2 k-1 � i=1 ∞ � j=i Γ j |� k-i | + k-1 � i=0 ∞ � j=i Γ j |� � k-i | + h(k),
while the loss of CBMPI satisfies

l k ≤ 2 k-2 � i=1 ∞ � j=i+m Γ j |� k-i-1 | + k-1 � i=0 ∞ � j=i Γ j |� � k-i | + h(k), where h(k) Δ = 2 � ∞ j=k Γ j |d 0 | or h(k) Δ = 2 � ∞ j=k Γ j |b 0 |.
Remark 1. A close look at the existing point-wise error bounds for AVI (Munos, 2007, Lemma 4.1) and API (Munos, 2003, Corollary 10) shows that they do not consider error in the greedy step (i.e., � � k = 0) and that they have the following form:

lim sup k→∞ l k ≤ 2 lim sup k→∞ k-1 � i=1 ∞ � j=i Γ j |� k-i |.
This indicates that the bound in Lemma 2 not only unifies the analysis of AVI and API, but it generalizes them to the case of error in the greedy step and to a finite horizon k. Moreover, our bound suggests that the way the errors are propagated in the whole family of algorithms VI/PI/MPI does not depend on m at the level of the abstraction suggested by Definition 1. 4

The next step is to show how the point-wise bound of Lemma 2 can turn to a bound in weighted L p -norm, which for any function f : S → R and any distribution µ on S is defined as �f � p,µ [START_REF] Munos | Error bounds for approximate policy iteration[END_REF][START_REF] Munos | Performance bounds in Lp-norm for approximate value iteration[END_REF]; [START_REF] Munos | Finite-time bounds for fitted value iteration[END_REF], and the recent work of [START_REF] Farahmand | Error propagation for approximate policy and value iteration[END_REF], which provides the most refined bounds for API and AVI, show how to do this process through quantities, called concentrability coefficients, that measure how a distribution over states may concentrate through the dynamics of the MDP. We now state a lemma that generalizes the analysis of [START_REF] Farahmand | Error propagation for approximate policy and value iteration[END_REF] to a larger class of concentrability coefficients. We will discuss the potential advantage of this new class in Remark 4. We will also show through the proofs of Thms. 1 and 3, how the result of Lemma 3 provides us with a flexible tool for turning point-wise bounds into L p -norm bounds. Thm. 3 in Appendix. D provides an alternative bound for the loss of AMPI, which in analogy with the results of [START_REF] Farahmand | Error propagation for approximate policy and value iteration[END_REF] shows that the last iterations have the highest impact on the loss (the influence exponentially decreases towards the initial iterations). Lemma 3 (Proof in Appendix. C). Let I and (J i ) i∈I be sets of positive integers, {I 1 , . . . , I n } be a partition of I, and f and (g i ) i∈I be functions satisfying

Δ = �� |f (x)| p µ(dx) � 1/p .
|f | ≤ � i∈I � j∈Ji Γ j |g i | = n � l=1 � i∈I l � j∈Ji Γ j |g i |.
Then for all p, q and q � such that 1 q + 1 q � = 1, and for all distributions ρ and µ, we have 4 Note however that the dependence on m will reappear if we make explicit what is hidden in the terms Γ j .

�f � p,ρ ≤ n � l=1 � C q (l) � 1/p sup i∈I l �g i � pq � ,µ � i∈I l � j∈Ji γ j ,
with the following concentrability coefficients

Cq(l) Δ = � i∈I l � j∈J i γ j cq(j) � i∈I l � j∈J i γ j ,
with the Radon-Nikodym derivative based quantity cq(j)

Δ = max π 1 ,••• ,π j � � � � d(ρPπ 1 Pπ 2 • • • Pπ j ) dµ � � � � q,µ . (14)
We now derive a L p -norm bound for the loss of the AMPI algorithm by applying Lemma 3 to the pointwise bound of Lemma 2.

Theorem 1 (Proof in Appendix. D). Let ρ and µ be distributions over states. Let p, q, and q � be such that 1 q + 1 q � = 1. After k iterations, the loss of AMPI satisfies

�l k �p,ρ ≤ 2(γ -γ k ) � C 1,k,0 q � 1 p (1 -γ) 2 sup 1≤j≤k-1 ��j� pq � ,µ (15)
+ (1 -γ k ) � C 0,k,0 q � 1 p (1 -γ) 2 sup 1≤j≤k �� � j � pq � ,µ + g(k),
while the loss of CBMPI satisfies

�l k �p,ρ ≤ 2γ m (γ -γ k-1 ) � C 2,k,m q � 1 p (1 -γ) 2 sup 1≤j≤k-2 ��j� pq � ,µ + (1 -γ k ) � C 1,k,0 q � 1 p (1 -γ) 2 sup 1≤j≤k �� � j � pq � ,µ + g(k), (16)
where for all q, l, k and d, the concentrability coefficients C l,k,d q are defined as

C l,k,d q Δ = (1 -γ) 2 γ l -γ k k-1 � i=l ∞ � j=i γ j cq(j + d),
with c q (j) given by Eq. 14, and g(k) is defined as

g(k) Δ = 2γ k 1-γ � C k,k+1 q � 1 p min � �d 0 � pq � ,µ , �b 0 � pq � ,µ � .
Remark 2. When p tends to infinity, the first bound of Thm. 1 reduces to

�l k �∞ ≤ 2(γ -γ k ) (1 -γ) 2 sup 1≤j≤k-1 ��j�∞ + 1 -γ k (1 -γ) 2 sup 1≤j≤k �� � j �∞ + 2γ k 1 -γ min(�d0�∞, �b0�∞). ( 17
)
When k goes to infinity, Eq. 17 gives us a generalization of the API (m = ∞) bound of Bertsekas & Tsitsiklis (1996, Prop. 6.2), i.e., lim sup

k→∞ �l k �∞ ≤ 2γ sup j ��j�∞ + sup j �� � j �∞ (1 -γ) 2 .
Moreover, since our point-wise analysis generalizes those of API and AVI (as noted in Remark 1), the L p -bound of Eq. 15 unifies and generalizes those for API [START_REF] Munos | Error bounds for approximate policy iteration[END_REF] and AVI [START_REF] Munos | Performance bounds in Lp-norm for approximate value iteration[END_REF].

Remark 3. [START_REF] Canbolat | Approximate) iterated successive approximations algorithm for sequential decision processes[END_REF] recently (and independently) developped an analysis of an approximate form of MPI. Also, as mentioned, the proof technique that we used is mainly based on that in [START_REF] Thiery | Performance bound for approximate optimistic policy iteration[END_REF]. While Canbolat & Rothblum (2012) only consider the error in the greedy step and [START_REF] Thiery | Performance bound for approximate optimistic policy iteration[END_REF] that in the value update, our work is more general in that we consider both sources of error -this is required for the analysis of CBMPI. [START_REF] Thiery | Performance bound for approximate optimistic policy iteration[END_REF] and [START_REF] Canbolat | Approximate) iterated successive approximations algorithm for sequential decision processes[END_REF] provide bounds when the errors are controlled in max-norm, while we consider the more general L p -norm. At a more technical level, Thm. 2 in [START_REF] Canbolat | Approximate) iterated successive approximations algorithm for sequential decision processes[END_REF] bounds the norm of the distance v *v k , while we bound the loss v *v π k .

If we derive a bound on the loss (using e.g., Thm. 1 in [START_REF] Canbolat | Approximate) iterated successive approximations algorithm for sequential decision processes[END_REF], this leads to a bound on the loss that is looser than ours. In particular, this does not allow to recover the standard bounds for AVI/API, as we managed to (c.f. Remark 2). Remark 4. We can balance the influence of the concentrability coefficients (the bigger the q, the higher the influence) and the difficulty of controlling the errors (the bigger the q � , the greater the difficulty in controlling the L pq � -norms) by tuning the parameters q and q � , given the condition that 1 q + 1 q � = 1. This potential leverage is an improvement over the existing bounds and concentrability results that only consider specific values of these two parameters: q = ∞ and q � = 1 in [START_REF] Munos | Performance bounds in Lp-norm for approximate value iteration[END_REF]; [START_REF] Munos | Finite-time bounds for fitted value iteration[END_REF], and q = q � = 2 in [START_REF] Farahmand | Error propagation for approximate policy and value iteration[END_REF]. Remark 5. For CBMPI, the parameter m controls the influence of the value function approximator, cancelling it out in the limit when m tends to infinity (see Eq. 16). Assuming a fixed budget of sample transitions, increasing m reduces the number of rollouts used by the classifier, and thus, worsens its quality; in such a situation, m allows to make a trade-off between the estimation error of the classifier and the overall value function approximation.

Finite-Sample Analysis of CBMPI

In this section, we focus on CBMPI and detail the possible form of the error terms that appear in the bound of Thm. 1. We select CBMPI among the proposed algorithms because its analysis is more general than the others as we need to bound both greedy and evaluation step errors (in some norm), and also because it displays an interesting influence of the parameter m (see Remark 5). We first provide a bound on the greedy step error. From the definition of � � k for CBMPI (Eq. 12) and the description of the greedy step in CBMPI, we can easily observe that ��

� k � 1,µ = L Π k-1 (µ; π k ).
Lemma 4 (Proof in Appendix. E). Let Π be a policy space with finite VC-dimension h = V C(Π) and µ be a distribution over the state space S. Let N be the number of states in D � k-1 drawn i.i.d. from µ, M be the number of rollouts per state-action pair used in the estimation of � Q k-1 , and π k = argmin π∈Π � L Π k-1 (� µ, π) be the policy computed at iteration k -1 of CBMPI. Then, for any δ > 0,

�� � k � 1,µ = L Π k-1 (µ; π k ) ≤ inf π∈Π L Π k-1 (µ; π) + 2(� � 1 + � � 2 ),
with probability at least 1δ, where

� � 1 (N, δ) = 16Qmax � 2 N � h log eN h + log 32 δ � , � � 2 (N, M, δ) = 8Qmax � 2 M N � h log eM N h + log 32 δ � .
We now consider the evaluation step error. The evaluation step at iteration k of CBMPI is a regression problem with the target (T

π k ) m v k-1 and a training set �� s (i) , � v k (s (i) )
�� n i=1 in which the states s (i) are i.i.d. samples from µ and � v k (s (i) ) are unbiased estimates of the target computed according to Eq. 7. Different function spaces F (linear or non-linear) may be used to approximate (T π k ) m v k-1 . Here we consider a linear architecture with parameters α ∈ R d and bounded (by L) basis functions {ϕ j } d j=1 , �ϕ j � ∞ ≤ L. We denote by φ :

X → R d , φ(•) = � ϕ 1 (•), . . . , ϕ d (•) � �
the feature vector, and by F the linear function space spanned by the features ϕ j , i.e., F = {f α (•) = φ(•) � α : α ∈ R d }. Now if we define v k as the truncation (by V max ) of the solution of the above linear regression problem, we may bound the evaluation step error using the following lemma.

Lemma 5 (Proof in Appendix. F). Consider the linear regression setting described above, then we have

�� k � 2,µ ≤ 4 inf f ∈F �(T π k ) m v k-1 -f � 2,µ + � 1 + � 2 ,
with probability at least 1δ, where

�1(n, δ) = 32Vmax � 2 n log � 27(12e 2 n) 2(d+1) δ � , �2(n, δ) = 24 � Vmax + �α * �2 • sup x �φ(x)�2 � � 2 n log 9 δ ,
and α * is such that f α * is the best approximation (w.r.t. µ) of the target function (T π k ) m v k-1 in F.

From Lemmas 4 and 5, we have bounds on �� � k � 1,µ and �� k � 1,µ ≤ �� k � 2,µ . By a union bound argument, we thus control the r.h.s. of Eq. 16 in L 1 -norm. In the context of Thm. 1, this means p = 1, q � = 1 and q = ∞, and we have the following bound for CBMPI: Theorem 2. Let d � = sup g∈F ,π � inf π∈Π L Π π � ,g (µ; π) and d m = sup g∈F ,π inf f ∈F �(T π ) m gf � 2,µ . With the notations of Thm. 1 and Lemmas 4 and 5, after k iterations, and with probability 1δ, the expected loss E µ [l k ] = �l k � 1,µ of CBMPI is bounded by

2γ m (γ -γ k-1 )C 2,k,m ∞ (1 -γ) 2 � dm + �1(n, δ 2k ) + �2(n, δ 2k ) � + (1 -γ k )C 1,k,0 ∞ (1 -γ) 2 � d � + � � 1 (N, δ 2k ) + � � 2 (N, M, δ 2k ) � + g(k).
Remark 6. This result leads to a quantitative version of Remark 5. Assume that we have a fixed budget for the actor and the critic B = nm = N M |A|m.

Then, up to constants and logarithmic factors, the bound has the form

�l k � 1,µ ≤ O � γ m � d m + � m B � + d � + � M |A|m B �
. It shows the trade-off in the tuning of m: a big m can make the influence of the overall (approximation and estimation) value error small, but that of the estimation error of the classifier bigger.

Summary and Extensions

In this paper, we studied a DP algorithm, called modified policy iteration (MPI), that despite its generality that contains the celebrated policy and value iteration methods, has not been thoroughly investigated in the literature. We proposed three approximate MPI (AMPI) algorithms that are extensions of the wellknown ADP algorithms: fitted-value iteration, fitted-Q iteration, and classification-based policy iteration.

We reported an error propagation analysis for AMPI that unifies those for approximate policy and value iteration. We also provided a finite-sample analysis for the classification-based implementation of AMPI (CBMPI), whose analysis is more general than the other presented AMPI methods. Our results indicate that the parameter of MPI allows us to control the balance of errors (in value function approximation and estimation of the greedy policy) in the final performance of CBMPI. Although AMPI generalizes the existing AVI and classification-based API algorithms, additional experimental work and careful theoretical analysis are required to obtain a better understanding of the behaviour of its different implementations and their relation to the competitive methods. Extension of CBMPI to problems with continuous action space is another interesting direction to pursue.

  and r

	(i) t and s	(i) t+1 are the reward
	and next state induced by this choice of action. From
	this rollout, we compute an unbiased estimate �

Note that the function space F and policy space Π are automatically defined by the choice of the regressor and classifier, respectively.

Here we used the same sampling distribution µ for both regressor and classifier, but in general different distributions may be used for these two components.

We may implement CBMPI more sample efficient by reusing the rollouts generated for the greedy step in the evaluation step.

Acknowledgments The second and third authors would like to thank French National Research Agency (ANR) under project LAMPADA n • ANR-09-EMER-007, European Community's Seventh Framework Pro-gramme (FP7/2007-2013) under grant agreement n • 231495, and PASCAL2 European Network of Excellence for supporting their research.


Panagiotis Papadakis
email: ppapadak@iit.demokritos.gr

Ioannis Pratikakis

Stavros Perantonis

Theoharis Theoharis

Efficient 3D shape matching and retrieval using a concrete radialized spherical projection representation

Keywords: 3D shape retrieval, Spherical harmonics, Descriptor's invariance

We present a 3D shape retrieval methodology based on the theory of spherical harmonics. Using properties of spherical harmonics, scaling and axial flipping invariance is achieved. Rotation normalization is performed by employing the continuous principal component analysis along with a novel approach which applies PCA on the face normals of the model. The 3D model is decomposed into a set of spherical functions which represents not only the intersections of the corresponding surface with rays emanating from the origin but also points in the direction of each ray which are closer to the origin than the furthest intersection point. The superior performance of the proposed methodology is demonstrated through a comparison against state-of-the-art approaches on standard databases.

Introduction

3D models have become an important part in modern computer graphics applications, such as computer-aided design, game development and film production while in several scientific domains, such as molecular biology, medicine, computational geometry and computer vision, a large part of scientific data is 3D data. The rapid evolution in graphics hardware and software development which has greatly facilitated 3D model acquisition, creation and manipulation, has given the opportunity to experience applications using 3D models not only to specialized users of the scientific community and the industrial domain, but also to common users. As the number of 3D models is continuously growing and thousands of models are already available from public and proprietary databases, ᭧ the problem in creating new 3D models may shift to the problem of searching for existing 3D models. Instead of creating 3D models from scratch, a more convenient and profitable approach is the use of existing 3D models. Thereupon, the development of efficient search mechanisms is needed for the retrieval of 3D models from large repositories.

3D model search and retrieval could be performed by using a textual description of the user's target which identifies the semantic meaning of the desired model or class of models. In this case, the user would explicitly describe the target, but such an approach is sensitive to the user's subjectivity factor which is not necessarily in agreement with the textual information which has been annotated to the target. Furthermore, this method is problematic as it requires individually annotating every model of a repository which is impractical due to the huge and continuously increasing number of existing 3D models. Therefore, content-based 3D shape retrieval methods are suited for search since they do not require any annotation while. They only require robust 3D shape feature extraction that can be applied automatically. In these methods, a shape descriptor is computed which represents the model and is consequently used at the matching stage. When 3D model comparison is performed, it is desired that shape descriptors are invariant under geometrical transformations such as translation, scaling, rotation and reflection. Thereafter, the discriminative power of these methods is highly affected by the invariance to these transformations, while extraction and comparison time also affect the performance, especially for real-time applications.

In this paper, we describe a 3D shape descriptor based on the theory of spherical harmonics. Translation invariance is achieved by using continuous principal component analysis (CPCA). The rotation invariance problem is alleviated by applying PCA on the face normals of the model (NPCA) and CPCA on the surface points of the model. The 3D model is decomposed into a set of spherical functions which represents not only the intersection points of the model's surface with rays emanating from the origin but also all points in the direction of each ray that are closer to the origin than the furthest intersection point. Spherical functions are then expressed by spherical harmonic coefficients by applying the spherical harmonics transform individually to each spherical function. Scaling and axial flipping invariance is achieved in the last stage of the shape descriptor's extraction, by using properties of spherical harmonic coefficients. The novelty of the proposed approach consists of the following: (i) Application of PCA for rotation normalization by using the face normals of the model and parallel use of CPCA on the surface points of the model; (ii) Scaling invariance is achieved by scaling the spherical functions proportionally to the model's size and by scaling the corresponding spherical harmonic coefficients to the unit L 1 norm; (iii) A new shape descriptor using spherical functions that represent not only the intersection points of the model's surface with rays emanating from the origin but also all points in the direction of each ray that are closer to the origin than the furthest intersection point.

The remaining of this paper is organized as follows: in Section 2, the state-of-the-art in the area of 3D model retrieval methods is presented. We summarize the approaches used to achieve rotation invariance in Section 3 and describe the proposed 3D shape descriptor. In Section 4, we present performance results and a comparison against start-of-the-art retrieval methods. Finally, in Section 5 conclusions and future work are discussed.

Related work

In this section we give an overview of existing start-of-theart 3D shape retrieval methods. We categorize these methods into the following: (i) spherical function based; (ii) 2D imageprojection based; (iii) volumetric features based; (iv) hierarchical structure based; (v) shape distributions based and (vi) hybrid methods.

In the first category, 3D space is parameterized using spherical functions. A popular example is the shape histograms descriptor proposed by Ankerst et al. [START_REF] Ankerst | Nearest neighbor classification in 3D protein databases[END_REF], which was developed primarily to classify proteins. In this method, 3D space can be divided into concentric shells, sectors, or both. Dividing 3D space into both sectors and shells gives the best results. For each part of the subdivision the model's shape distribution is computed, giving a sum of histogram bins which comprises the final shape descriptor. In [START_REF] Körtgen | 3D shape matching with 3D shape contexts[END_REF] Koertgen et al. make use of 3D shape contexts which are histograms each one identifying a surface point by the relative coordinates of the remaining surface points. Vranic et al. [START_REF] Vranic | 3D model retrieval[END_REF] proposed the ray-based descriptor which characterizes a 3D model by a spherical extent function capturing the furthest intersection points of the model's surface with rays emanating from the origin. Spherical harmonics and moments representation of the spherical extent function were tested in [START_REF] Saupe | 3d model retrieval with spherical harmonics and moments[END_REF], indicating better performance in the first case. In the same context, a more general approach was introduced by Vranic et al. [START_REF] Vranic | 3D model retrieval[END_REF][START_REF] Vranic | An improvement of rotation invariant 3D shape descriptor based on functions on concentric spheres[END_REF] known as the layered depth spheres descriptor (also known as radialized spherical extent function). Here, the 3D model is described by a spherical function which decomposes the model into a sum of concentric shells and gives the maximal distance of the model from the center of mass as a function of angle and the radius of the equivalent shell. The spherical function is represented by spherical harmonics coefficients. The extended Gaussian image (EGI) proposed by Horn et al. [START_REF] Horn | Extended Gaussian images[END_REF] has also been used as a shape representation method, where a spherical function correlates the faces of a 3D model with values on a sphere of radius 1 (Gaussian sphere). Specifically, the Gaussian sphere is discretized into a finite set of points. Then, each point is assigned a value proportional to the surface area of the model having the same direction as the unit vector defined by the specific point. A variation of this basic idea is the complex EGI proposed by Kang et al. [START_REF] Kang | Determining 3-D object pose using the complex extended Gaussian image[END_REF] where each EGI cell is described by a complex number which captures both rotation and translation information. Novotni et al. [START_REF] Novotni | 3D Zernike descriptors for content based shape retrieval[END_REF] presented the 3D Zernike descriptor, generalizing the 2D Zernike descriptor. The descriptor is an extension of spherical harmonics descriptors because Zernike functions are spherical harmonic functions modulated by appropriate radial functions. Norms of vectors consisting of Zernike moments are used to comprise the shape descriptor in order to achieve rotation invariance. In [START_REF] Daras | Efficient 3-D model search and retrieval using generalized 3-D radon transforms[END_REF], Daras et al. introduces the generalized radon transform, where the descriptor is extracted by combining the radial integration transform which captures the radial information of a 3D model and the spherical integration transform which captures its spherical information.

In the second category of 3D shape retrieval methods, 2D images-projections of a 3D model are used for feature extraction and comparison. Chen et al. [START_REF] Chen | On visual similarity based 3D model retrieval[END_REF] proposed the light field descriptor, which is presented in the literature as the method with the best discriminative power so far, but having long extraction and comparison times as a disadvantage. The shape descriptor is generated by projecting the 3D model uniformly from the vertices of a dodecahedron and the distance between two models is computed by the minimal L 1 distance between all pairs of vertices under all possible rotations of the two dodecahedrons. Each projection is composed of a combined feature vector of Zernike moments and Fourier coefficients. Another method where 2D shape matching can support 3D shape matching is the depth-buffer descriptor presented in [START_REF] Vranic | 3D model retrieval[END_REF] by Vranic, where six projections of the model are created, one for each side of a cube which encloses the model, storing the maximum and minimum depth values of the model at each side.

The six depth buffers can then be used to produce shape descriptors in the spatial domain or in the spectral domain. In the same work, the silhouette-based descriptor is proposed, which is produced by projecting the 3D model to the Cartesian planes, leading to three silhouettes. The contours of the silhouettes are computed by two alternative approaches and are used accordingly at the matching stage. Ohbuchi et al. [START_REF] Ohbuchi | Retrieving 3D shapes based on their appearance[END_REF] proposed the multiple orientation depth Fourier transform descriptor. In this method, the model is projected from 42 viewpoints to cover all possible view aspects of a 3D model. The respective depth buffers are then transformed to the r-coordinate system and a Fourier transform is applied, producing a set of coefficients for each depth buffer. To compare two models, all possible pairs of coefficients are compared which inevitably increases comparison time. Zarpalas et al. [START_REF] Zarpalas | 3D model search and retrieval using the spherical trace transform[END_REF] introduced a 3D shape descriptor called the spherical trace transform, which is an extension to 3D shape matching of the 2D trace transform. Among a variety of 2D features which are computed over a set of planes intersecting the volume of a 3D model, one of the components of the shape descriptor are 2D Krawtchouk moments which improve significantly the overall performance of the proposed method.

In the case of methods that belong to the third category, initially the 3D model is preprocessed to obtain a volumetric representation which is then used at the feature extraction stage. This is done by dividing 3D space into volume elements (voxels), and then assigning to each voxel a value according to a function which considers the model's relative position. Many variations exist concerning the method used to decide what value is to be assigned to each voxel. The Gaussian (exponentially decaying) Euclidean distance transform proposed by Kazhdan et al. [START_REF] Kazhdan | Rotation invariant spherical harmonic representation of 3D shape descriptors[END_REF] (GEDT) is one of best performing variations. The GEDT descriptor can be expressed by norms of spherical harmonic frequencies which give a rotation invariant representation, reducing the shape descriptor's storage cost and comparison time. The reflective symmetry descriptor proposed by Kazhan et al. [START_REF] Kazhdan | A reflective symmetry descriptor[END_REF], computes a voxelized representation of a 3D model which is then used to compute the reflective symmetries of the model with respect to all planes passing through its center of mass. Suzuki et al. [START_REF] Suzuki | A similarity retrieval of 3D polygonal models using rotation invariant shape descriptors[END_REF], introduced a shape descriptor based on equivalence classes. After enclosing the model in a 3D grid, cells which follow the same paths under 90 • rotations are organized into sets called equivalence classes. Each of these classes is represented by a value giving the density of the initial point cloud inside the cells of the given class. Tangelder et al. [START_REF] Tangelder | Polyhedral model retrieval using weighted point sets[END_REF], make use of weighted point sets. The set of points is selected by a 3D grid which encloses the model and three methods are proposed to select a point at each non-empty grid cell and assign a weight to the specific point. For the comparison of two shape signatures a variation of the earth mover's distance is adopted.

The fourth category consists of methods which use hierarchical 3D model representations based on skeleton graphs. Hilaga et al. [START_REF] Hilaga | Topology matching for fully automatic similarity estimation of 3D shapes[END_REF] introduced the multiresolution Reeb graph, which represents a 3D model's topology and skeletal structure at various levels of detail. In [START_REF] Zhang | Retrieving articulated 3-D models using medial surfaces and their graph spectra[END_REF] Zhang et al. consider the use of medial surfaces to compute an equivalent directed acyclic graph of a voxelized model. Sundar et al. [START_REF] Sundar | Skeleton based shape matching and retrieval[END_REF] also proposed skeleton based shape matching. They first obtain the volume of a 3D model which passes through a thinning process producing a set of skeletal points, which finally form a directed acyclic graph by applying the minimum spanning tree algorithm. Cornea et al. [START_REF] Cornea | 3D object retrieval using many-to-many matching of curve skeletons[END_REF], propose the use of curve skeletons produced by the application of the generalized distance field to the volumetric representation of a 3D model. Matching between two skeletons is performed by using the earth mover's distance. Generally, graph-based methods for 3D model matching are mostly useful for matching articulated models but are problematic in the use of efficient similarity measures.

Methods based on statistical measures belong to the last category. The shape distributions proposed by Osada et al. [START_REF] Osada | Matching 3D models with shape distributions[END_REF][START_REF] Osada | Shape distributions[END_REF] is a shape descriptor where several shape characteristics can be measured for a random set of points belonging to the model, according to the selection of an appropriate shape function. There are variations of this method according to the characteristic that is measured, e.g. in the D2 shape distribution descriptor the distance between two random surface points is computed. Ohbuchi et al. [START_REF] Ohbuchi | Shape-similarity search of 3D models by using enhanced shape functions[END_REF], proposed enhanced shape functions, namely the angle-distance histogram or the absolute angle-distance histogram for inconsistently oriented meshes, which are extensions of the D2 shape distribution descriptor and have better performance. In [START_REF] Ohbuchi | Shape-similarity comparison of 3D models using alpha shapes[END_REF], the same author proposes a combination of the absolute angle-distance function with the alpha shape representation of 3D models. Zaharia et al. [START_REF] Zaharia | 3D shape-based retrieval within the MPEG-7 framework[END_REF] presented the 3D shape spectrum descriptor which is within the MPEG-7 framework. This shape representation is the histogram of the shape index values calculated for the entire mesh of a 3D model, where each shape index value is the angular representation of the vector composed of the first and second principal curvature at a specific surface point.

Finally, the hybrid descriptor proposed by Vranic [START_REF] Vranic | 3D model retrieval[END_REF] is a combination of three different methods. The depth buffer, silhouette and ray-based descriptors described earlier, are used together into one descriptor trying to capture as much discriminating information for a 3D model as possible.

The proposed methodology

The overall scheme

In this section, we provide an overview of the proposed methodology that will be detailed in Sections 3.2-3.4. The main steps for the 3D shape descriptor's extraction procedure as well as the matching scheme are depicted in Fig. 1.

Initially, the model is translated so that its center of mass coincides with the origin, as shown in Fig. 1(a). Then we apply in parallel two different alignment methods, namely CPCA and the PCA on the normals of the model's faces (NPCA), to alleviate the problem of rotation invariance. This gives two aligned versions of the translated model (Fig. 1(b)) which are hereafter processed in the same way and will finally give two 3D shape descriptors. In the next stage, each alignment variant is expressed by a set of spherical functions f 1 , f 2 , . . . , f N with increasing radius giving the intersections of the model's surface with rays emanating from the origin (Fig. 1(c)). This is equivalent to producing spherical projections of approximately equidistant parts of the model at different radii. In the sequel, the set of spherical functions is processed to find the furthest intersection points from the origin on each ray. If the model is viewed in the direction of a specific ray, then the furthest intersection point obscures the part of the model along the ray which is closer to the origin. We may assume that the obscured part belongs to the 3D model because we perceive a 3D model as a concrete entity and not as a 3D surface composed of a set of faces. This may not be true for concave objects but our experiments showed that this does not influence the discriminative power of the shape descriptor. Thereafter, these parts are considered to belong to the model giving a new 3D model representation (Fig. 1(d)). For each spherical function, the spherical harmonics transform is applied producing a set of complex coefficients (Fig. 1(e)). In the last stage, spherical harmonic coefficients are scaled to the unit L 1 norm (Fig. 1(f)) which guarantees scaling invariance and properties of spherical harmonics are used to achieve axial flipping invariance. From here on, we will call the proposed descriptor the concrete radialized spherical projection (CRSP) descriptor.

Each 3D model is represented by the concatenation two shape descriptors corresponding to the aligned versions using CPCA and NPCA and of a set of complex spherical harmonic coefficients. To compare two shape descriptors, the L 1 distance between the corresponding parts is computed and the minimum distance is taken as the similarity measure between the two models. The notion of taking the minimum distance is based on the expectation that the best establishment of correspondences between two models is achieved, when the difference between the corresponding shape descriptors is minimum.

Translation and rotation normalization

In this section, we detail how we address the problem of translation and rotation invariance. Prior to the translation and rotation normalization procedures, a preprocessing step is applied to obtain a representation of a 3D model, consisting only of triangles. A double triangle check is also performed to find identical triangles within the polygonal model. These triangles are recognized as having the same vertices but appearing in different order within a 3D model file (e.g. the triangles with vertex indices (10, 422, 87) and (422, 10, 87) are identical). We found a large number of such triangles within the 3D model datasets tested in Section 4. These triangles cannot be identified when a 3D model is rendered, but if they are taken into account they may affect the model's translation, rotation and scale normalization.

The model is translated to its center of mass m which is computed using CPCA [START_REF] Vranic | 3D model retrieval[END_REF][START_REF] Vranic | Tools for 3D-object retrieval: Karhunen-Loeve transform and spherical harmonics[END_REF].

If v = [v x , v y , v z ] T is a surface point, then using barycentric coordinates, v is given by v = a • A i + b • B i + (1 -a -b) • C i , (1)
where A i , B i and C i are the vertices of triangle i, i = 1, 2, . . . , N, where N is the number of triangles of the model and a, b are the barycentric coordinates. Thus m is equal to:

m = 1 E v∈T i v ds = 1 E N i=1 E i • (A i + B i + C i ) 3 , ( 2
)
where E is the total surface area of the model and E i is the surface area of triangle T i . The detailed derivation of Eq. ( 2) is given in Appendix A. When a shape descriptor is based on the establishment of correspondences between models, the rotation invariance property is of great importance. Translation and scaling invariance are also important and a shape descriptor can be invariant under these transformations using efficient methods. However, to achieve rotation invariance using a single method for a wide range of model classes, has proven to be a very complicated problem. This is because it is not obvious which shape characteristic we should choose and the way it should be used when normalizing for rotation, so that it is suitable for every class. Several approaches have been used belonging to one of the following categories: (i) alignment and (ii) rotation invariance by definition.

In category (i), each model is normalized for rotation by rotating the model into a canonical coordinate frame with the expectation that when comparing models of the same class these models will be similarly aligned. However, as mentioned previously, 3D models are spread over a wide variety of classes, making it very difficult to find an alignment method for models of every possible class. The prominent alignment method is PCA, also known as the Karhunen-Loeve transform [START_REF] Theodoridis | Pattern Recognition[END_REF]. The PCA algorithm generates mutually uncorrelated features taking a set of vectors as input. The set of vertices of the polygonal model [START_REF] Vranic | 3D model retrieval[END_REF] or the centroids of the model's triangles [START_REF] Paquet | Description of shape information for 2-D and 3-D objects[END_REF] have been used as input to the PCA. In [START_REF] Vranic | 3D model retrieval[END_REF][START_REF] Vranic | Tools for 3D-object retrieval: Karhunen-Loeve transform and spherical harmonics[END_REF] the so-called CPCA is presented. Another approach to normalize a 3D model for rotation (similar to PCA), is the use of singular value decomposition (SVD) [START_REF] Theodoridis | Pattern Recognition[END_REF]. In [START_REF] Osada | Matching 3D models with shape distributions[END_REF][START_REF] Elad | Contend based retrieval of VRML objects-an iterative and interactive approach[END_REF], the SVD of the covariance matrix of the 3D model is computed and the unitary matrix is applied to the model for rotation normalization.

In the category (ii), rotation invariance is achieved by the definition of the shape descriptor. These descriptors are invariant under rotation, but they may discard information regarding the 3D model which can be discriminative. Descriptors based on spherical harmonics, Zernike moments and shell histograms are examples of representation methods able to achieve rotation invariance by definition.

Instead of just using one method to achieve rotation invariance, combinations of different approaches may be considered. This may improve the performance of a shape descriptor at the cost of increased complexity.

To alleviate the rotation invariance problem, we apply two alternative alignment procedures in parallel. In the first procedure, we use CPCA to compute the covariance matrix C of the model's surface, as in the following:

C = 1 E v∈M (v -m)(v -m) T ds = 1 E N i=1 v∈T i (v -m)(v -m) T ds ⇒ C = 1 12E N i=1 E i • [f(A i ) + f(B i ) + f(C i ) + 9 • f((A i + B i + C i )/3)], ( 3
)
where M is the set of vertices of the polygonal model and

f(v) = (v -m)(v -m) T .
A detailed derivation of Eq. ( 3) is given in Appendix A; this is not included in [START_REF] Vranic | 3D model retrieval[END_REF][START_REF] Vranic | Tools for 3D-object retrieval: Karhunen-Loeve transform and spherical harmonics[END_REF] where CPCA was introduced.

The eigenvectors of C are then computed, which represent the principal directions of the model's surface area and form an orthonormal basis for R 3 space. These vectors are then sorted according to the order of descending (or ascending) eigenvalues and are finally normalized to the Euclidean unit norm. A rotation matrix R is formed having as columns the scaled eigenvectors in decreasing (or ascending) order and R is applied to the model so that the principal directions of the model coincide with the coordinate axes. After the application of R, the model is rotated to a canonical position as shown in Fig. 2.

In the second alignment procedure we compute the covariance matrix by using NPCA. If n i is the normal vector of triangle T i and m Nl is the average face normal, then the covariance matrix is computed as

C = 1 2E N i=1 (E i (n i -m Nl ) • (n i -m Nl ) T + E i (-n i -m Nl ) • (-n i -m Nl ) T ). ( 4
)
Each triangle is counted twice, using its normal and opposite to the normal direction. This is done to avoid possible variations in vertex ordering at the model triangles, that affects the direction's sign and because we may assume that every triangle has two sides with opposite directions. Adopting this, we observe that the average face normal m Nl becomes

m Nl = 1 2E N i=1 (E i • n i + E i • (-n i )) = 0. (5)
The above equation indicates that m Nl is always the zero vector, thus C is now given by After the computation of C, the PCA algorithm is applied to find and align the principal axes of the model with the coordinate axes. Our experiments indicate that there are models where NPCA gives a more consistent alignment than CPCA and vice versa. To justify this, we give some examples in Fig. 3. Motivated by this observation, we use both methods, expecting to achieve better overall results. This is justified by the experimental results given in Section 4. After completion of this stage we have two alignment variants of the initial model which will finally lead to two 3D shape descriptors.

C = 1 2E N i=1 (E i • n i • n T i + E i • (-n i ) • (-n T i )) = 1 E N i=1 E i • n i n T i . (6)

Representation of a polygonal model as a set of spherical functions

After alignment, the model's surface is represented by a set of spherical functions.

A spherical function describes a bounded area of the model, defined by a lower and an upper radius which delimit a spherical shell. This shell is the volume in which the underlying surface of the model is represented by the equivalent spherical function points. The region corresponding to the rth spherical function is the spherical shell defined in the interval [l r , l r+1 ) where the l r bound is given by

l r = 0, r = 1, (r -0.5) • M • d avg /N, 1 < r N + 1, ( 7
)
where M is a weighting factor which determines the radius of the largest sphere (we set M = 3), N is the number of spherical functions (we set N =24) and d avg is the average distance of the model's surface from the center of mass. Multiplying by M = 3 ensures that the model is adequately enclosed by the sphere with the largest radius and only outlying parts of the model may be excluded from the furthest boundary. The average distance is computed as

d avg = C 11 + C 22 + C 33 , (8)
where C ij is the (i, j )th element of the covariance matrix, com- We compute the intersections of the triangles of the aligned model with rays emanating from the origin in the ( j , k ) directions denoted as ray( j , k ), where j = (2j + 1) /(4B), k = 2 k/(2B), j , k = 0, 1, . . . , 2B -1 and B is the sampling bandwidth. Inters(i, j , k ) denotes the distance from the origin of the ith intersection of ray( j , k ) with the model's surface. If ray( j , k ) does not intersect the model's surface then Inters(0, j , k ) = 0.

Let mxin(r, j , k ) be the distance of the furthest intersection point along ray( j , k ), within the boundaries of the rth spherical function, given by mxin(r, j , k ) = max{Inters(i, j , k )|l r Inters(i, j , k ) < l r+1 , i = 0, 1, . . .}.

The 3D model is represented as a set of N spherical functions f r : S 2 → [l r , l r+1 ) ∪ {0}, r = 1, 2, . . . , N, where S 2 denotes the points ( j , k ) on the unit sphere and f r is given by

f r ( j , k ) = mxin(r, j , k ). ( 10
)
The decomposition of the 3D model into the spherical functions f r ( , ) is similar to projecting approximately equidistant parts of the 3D model to concentric spheres of increasing radius. The set of spherical functions representing the 3D model as given by Eq. ( 10) is depicted in Fig. 5, where Fig. 6. Intersections of a model's surface with rays emanating from the origin (circled dots) and points closer to the origin than the furthest intersection point at each ray (single dots).

for each spherical function point f r ( j , k ) = 0, a black dot is used to represent the equivalent intersection.

In the following, the spherical functions are expanded to include additional information. Let mxin(D, j , k ) = 0 (corresponding to the furthest intersection point on ray( j , k )), be assigned to the spherical function point f D ( j , k ). Then all f r ( j , k ), for r < D are to:

f r ( j , k ) = r • M • d avg /N ∀r < D. ( 11
)
Extra points are thus included in the spherical functions as depicted in Fig. 6. Practically, if the model is viewed from the ( j , k ) direction then the part of the model along ray( j , k ) which is closer to the origin than mxin(D, j , k ) is occluded by the intersection point corresponding to mxin(D, j , k ). We may assume that the occluded part actually belongs to the model because we perceive 3D models as concrete entities. Adopting this assumption, a more concrete representation of the 3D model is produced which highly increases discriminative power as shown in Section 4. Therefore, taking into account Eqs. ( 10), [START_REF] Chen | On visual similarity based 3D model retrieval[END_REF], the new spherical functions f r are

f r ( j , k ) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ mxin(r, j , k ), if mxin(D, j , k ) = 0 and r = D, r • M • d avg /N, if mxin(D, j , k ) = 0 and r < D, 0, otherwise. (12)
These spherical functions give a model representation as shown in Fig. 7, where black dots are used to illustrate spherical function points f r ( j , k ) = 0. In the sequel, we expand the N spherical functions to their spherical harmonic representation. A detailed description of the spherical harmonics' formulation is given in Appendix B. Thus, having formed N spherical functions, the spherical harmonics transform (SHT) (Appendix B, Theorem 1) is applied for each function. For the computation of SHT, some researchers have used the SpharmonicKit By applying the SHT, B 2 complex coefficients are produced for each spherical function. Higher order coefficients represent finer details of the model as well as possible noise. Instead of taking into account all B 2 coefficients, we used the first L bands of coefficients (up to the lth degree, we set L bands = 16). For all N spherical functions there are N * L 2 bands coefficients, where L bands B. Using Theorem 1 (Appendix B) and the fact that f r ( j , k ) is a real-valued function, the following property holds between the coefficients of positive and negative order:

f (l, m) = (-1) m f (l, -m), (13)
where x denotes the complex conjugate of x. This means that only the positive order coefficients, which are N * (L bands * (L bands + 1))/2 in total, are necessary to find the L 1 distance between two shape descriptors as will be proven in Section 3.4. We next show that the shape descriptor is invariant to possible axial flipping using the spherical harmonics' properties below:

(i) if the model is reflected with respect to the x-axis (yz plane) then: f (l, m) ⇒ f (l, m); (ii) if the model is reflected with respect to the y-axis (xz plane) then: f (l, m) ⇒ (-1) l+m f (l, m); (iii) if the model is reflected with respect to the z-axis

(xy plane) then: f (l, m) ⇒ (-1) m f (l, m).
These properties indicate that only the sign of the real and imaginary part of a coefficient may change, so we store the absolute values of the imaginary parts for all l, m and the absolute values of the real parts when l or m is odd. If l and m are both even, then the sign of the real part does not change and we store the actual values.

Finally, the spherical harmonic coefficients are scaled to the unit L 1 norm. From the definition of spherical harmonics, spherical harmonic coefficients are exactly proportional to the model's scale. So scaling invariance is achieved by scaling the coefficients which is done in constant time because we sample all models with the same bandwidth.

Similarity measure

Each spherical harmonic coefficient is a complex number. To compare two shape descriptors, the L 1 distance is computed.

Let f (l, m) = Re + j Im, f (l, m) = Re + j Im be two corresponding coefficients of two different shape descriptors. The L 1 distance between them is

| f (l, m) -f (l, m)| = (Re -Re ) 2 + (Im -Im ) 2 . ( 14
)
Eq. ( 14) is computed for all coefficients where l L bands and m 0, due to the property derived from Eq. ( 13) that leads to the following:

| f (l, m) -f (l, m)| = |(-1) m f (l, -m) -(-1) m f (l, -m)| = | f (l, -m) -f (l, -m)|. ( 15
)
Thus only the positive order coefficients are necessary. The descriptor consists of two sets of coefficients corresponding to the two aligned versions of the model (using CPCA and NPCA). The comparison between two models is done between the corresponding sets. Consequently, the CPCA aligned model is compared with the CPCA aligned version of the second model, and similarly for the NPCA version. The comparison giving the minimum distance sets the final difference. The notion of taking the minimum distance is based on the expectation that the best establishment of correspondences between two models is achieved, when the difference between the shape descriptors is minimum.

Experimental results

In this section, we give performance results for the proposed shape descriptor and compare it against state-of-the-art methods. The evaluation was run in the following 3D model databases: The Princeton Shape Benchmark (PSB) (both training and test set) [32], the classified models of the National Taiwan University database (NTU) [33], the MPEG-7 dataset [34], and the classified models of the CCCC dataset [START_REF]Dejan Vranic's 3D Search Engine[END_REF]. From the NTU and CCCC datasets only the classified models were used because the entire datasets contain a large number of unclassified models thus, using the whole datasets would not give accurate results. Table 1 shows the number of categories and the total number of models for each dataset.

We compare the proposed descriptor with the following shape retrieval methods which are considered state-of-the-art with respect to their discriminative power: (i) hybrid (DSR472) [START_REF] Vranic | 3D model retrieval[END_REF]; (ii) depth buffer-based (DB) [START_REF] Vranic | 3D model retrieval[END_REF]; (iii) silhouette-based (SIL) [START_REF] Vranic | 3D model retrieval[END_REF]; (iv) the ray-based using spherical harmonics (RSH) [START_REF] Vranic | Tools for 3D-object retrieval: Karhunen-Loeve transform and spherical harmonics[END_REF]; (v) light-field (LF) [START_REF] Chen | On visual similarity based 3D model retrieval[END_REF]; (vi) the spherical harmonic representation of the Gaussian Euclidean distance transform descriptor (SH-GEDT) [START_REF] Kazhdan | Rotation invariant spherical harmonic representation of 3D shape descriptors[END_REF].

Descriptors (i)-(iv) were tested using tools provided in [START_REF]Dejan Vranic's 3D Search Engine[END_REF], the LF descriptor was produced by using the executable provided at [33] and the SH-GEDT descriptor was produced by using the executable provided at [36].

Our quantitative evaluation is based on the following measures:

• Precision-Recall plots (P-R): For a query model which belongs to a class C, recall is the percentage of models of class C that are retrieved and precision is the proportion of models that belong to class C to the number of retrieved models. • Nearest neighbor: Given a query model, its nearest neighbor is its closest match. The nearest neighbor measure is the percentage of queries where the closest match belongs to the query model's class. • First tier, second tier: Given a query model belonging to class C with |C| number of models, the first tier measure is the proportion of the first |C| -1 retrieved models that belong to class C. The second tier measure is the proportion of the first (2|C| -1) retrieved models that belong to class C. We average these values for all models in each dataset. • Discounted cumulative gain (DCG): For a query model belonging to class C, a list G is produced where G i is equal to 1 if the ith match belongs to class C, otherwise G i = 0. Then DCG i is given by

DCG i = G 1 , i = 1, DCG i-1 + G i log 2 i otherwise (16)
and the overall score is given by

DCG = DCG k 1 + |C| i=2 log 2 i , ( 17
)
where k is the number of models of the dataset. For all the above quantitative measures, the best score is 100%. Fig. 8. P-R plot for the Princeton Shape Benchmark test dataset. A significant gain in discriminative power is achieved when NPCA is used along with CPCA to alleviate the rotation invariance problem and when the concrete representation is adopted (Eq. ( 11)). In the P-R plot of Fig. 8, we show the significant gain in discriminative power when Eq. ( 11) is adopted to obtain a more concrete representation of a 3D model and when CPCA and NPCA are used together to alleviate the rotation invariance problem.

In Fig. 9 the P-R plot comparing all descriptors on all datasets is given. Precision changes along the vertical axis and recall along the horizontal axis. For each method, p 50 and p 100 values are given in the legend of each figure, indicating average precision for recall 5-50% and 5-100%, respectively.

Comparing the plots as well as the p 50 and p 100 values, we conclude that the CRSP descriptor has an overall superior performance. The difference is greater for recall up to 50%, which is the portion of the models of a class that a user is mostly interested in. In addition, the gradient of the precision curve does not seem to change significantly for recall higher than 50%.

In Table 2, we give the corresponding scores of the proposed descriptor for the PSB test dataset and compare it to the LF and SH-GEDT descriptors. Their scores were taken from the Princeton Shape Benchmark [START_REF] Shilane | The Princeton Shape Benchmark[END_REF] and we used the tools provided in [32] to produce the scores of our descriptor. The top performance of the CRSP descriptor is thus justified by the quantitative measures as well.

We believe that a difference over 2% is not trivial when top performing methods are compared, thus the CRSP de- scriptor has a significant gain over the LF and the SH-GEDT descriptors.

In Table 3 the scores of the CRSP descriptor are given for the remaining datasets.

In Table 4, we give the quantitative measure scores of CRSP for each class of the "coarse2Test" and "coarse2Train" classifications of the Princeton Shape Benchmark and the respective DCG scores for the SH-GEDT and LF descriptors. The cardinality of each class is given in parenthesis next to the class name. It is evident that while the performance of CRSP may vary among different classes, it achieves a better overall performance than LF and SH-GEDT. As it is shown in Table 4, the best performance is achieved for the class "vehicle" having the highest cardinality, while the worst performance is for the "building" class having the lowest cardinality. The higher cardinality of a class the more significant the results are.

In Table 5, we have split the "plant" class of the coares2Test classification into its five subclasses according to the test classification, to measure how performance is affected in the case of different degrees of intra-class variability. All other classes (vehicle, animal, household, building and furniture) in the coarse2Test classification keep the same granularity thus the performance of the corresponding methods in these classes remains the same.

The CRSP descriptor has the additional advantage of having very short extraction and comparison time. Using a computer with an AMD processor at 1.4 GHz, the average extraction time (including searching and removing identical triangles) was less than one second and the average comparison time was 0.16 ms. The average number of triangles per model was 8180 and the average number of vertices was 4427. We believe that the short extraction and comparison time is a significant advantage of our descriptor that enables very efficient online search and retrieval of 3D models. Furthermore, they enable extensions and future improvement (such as combinations with other descriptors).

In Fig. 10, we show examples of retrieval results for different queries within the PSB test dataset. For each query (top left model), the first 16 matches are shown using the CRSP and the LF descriptors.

Comparing the retrieved models for each query, we can conclude that the quality of the results depends on the cardinality of the query model's class, the intra-class variation as well as the inter-class variation. The third factor is particularly important as there may be classes closely related to each other. In such a setting, mixed retrieved results from these classes will have little impact on the user's dissatisfaction, while the opposite is true when mixed retrieved models belong to classes with large inter-class variation. The above query results are from the Princeton Shape Benchmark 3D model dataset where models have been classified by the designers according to similarity in the shape's functionality. The CRSP descriptor is based on the establishment of correspondences between 3D models and measures the similarity between models by taking into account the shape's geometry. Thus, the quality of the results is even better when models are classified according to this criterion.

Conclusion

We propose a new methodology for content based 3D shape retrieval. It uses CPCA along with PCA applied on the model's face normals in order to deal with the rotation invariance problem. The model is decomposed into a set of spherical functions representing intersections of the model's surface with rays emanating from the origin. Parts closer to the origin than the furthest intersection point at each ray are included in the spherical functions, which improves the descriptor's discriminative power. Spherical functions are represented by spherical harmonics and their properties are used to achieve scaling invariance and axial flipping invariance.

We extensively evaluated the proposed 3D shape descriptor on various 3D model datasets and compared it against state-ofthe-art methods. The evaluation results lead to the conclusion that the CRSP descriptor has very high discriminative power and is ranked first among the other methods that we compared In future work we aim to further improve the quality of the results by considering the user's relevance feedback. In that framework, the user would mark relevant models or assign a measure of relevance to some of the results. Using this information the similarity between the query and the models of the database would be re-evaluated and a more accurate list of results would be returned.

v = a • A i + b • B i + (1 -a -b) • C i ,
where A i , B i and C i are the vertices of triangle T i , i =1, 2, . . . , N where N is the number of triangles of the model and a, b are the barycentric coordinates. If E is the total surface area of the 3D model and E i is the surface area of triangle T i , then using CPCA the model's center of mass m is equal to =

m = 1 E v∈T i v ds = 1 E N i=1 (aA i + bB i + (1 -a -b)C i ) ds = 2 E N i=1 E i 1 0 1-a 0 (aA i + bB i + (1 -a -b)C i ) db da = 2 E N i=1 E i 1 0 abA i + b 2 2 B i + b -ab - b 2 2 C i 1-a 0 da
E N i=1 E i 1 0 a(1 -a)A i + (1 -a) 2 2 B i + (1 -a) -a(1 -a) - (1 -a) 2 2 C i da = 2 E N i=1 E i 1 0 aA i -a 2 A i + 1 2 B i -aB i + a 2 2 B i + 1 2 C i -aC i + a 2 2 C i da = 2 E N i=1 E i a 2 2 (A i -B i -C i ) - a 3 3 A i - 1 2 B i - 1 2 C i + 1 2 a(B i + C i ) 1 0 = 2 E N i=1 E i 1 2 (A i -B i -C i ) - 1 3 A i - 1 2 B i - 1 2 C i + 1 2 (B i + C i ) ⇒ m = 1 E N i=1 E i • (A i + B i + C i ) 3 . 2
In

•(aA i + bB i + (1 -a -b)C i -m) T db da = 2 E N i=1 E i 1 0 ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 1-a 0 ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ a 2 A i A T i + b 2 B i B T i + (1 -2a -2b + 2ab + a 2 + b 2 )C i C T i + ab(A i B T i + B i A T i ) + (a -a 2 -ab)(A i C T i + C i A T i )+ (b -ab -b 2 )(B i C T i + C i B T i ) -aA i m T -bB i m T -(1 -a -b)C i m T - amA T i -bmB T i -(1 -a -b)mC T i + mm T db ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ da = 2 E N i=1 E i 1 0 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ (a 2 -a 3 )A i A T i +
(1a) 3 3

B i B T i + C i C T i + B i C T i 2 + C i B T i 2 +a (1 -a) 2 2 (A i B T i + B i A T i + A i C T i + C i A T i ) -a(1 -a)(A i m T + mA T i ) - (1 -a) 2 2 (B i m T + C i m T + mB T i + mC T i ) + (1 -a)mm T ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ da = 1 12E N i=1 E i 2A i A T i + 2B i B T i + 2C i C T i + B i C T i + C i B T i + A i B T i + B i A T i + A i C T i + C i A T i -4A i m T -4mA T i -4B i m T -4C i m T -4mB T i -4mC T i + 12mm T = 1 12E N i=1 E i ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ A i A T i -A i m T -mA T i + mm T + B i B T i -B i m T -mB T i + mm T +C i C T i -C i m T -mC T i + mm T + A i A T i + A i B T i + A i C T i -A i 3m T +B i A T i + B i B T i + B i C T i -B i 3m T + C i A T i + C i B T i + C i C T i -C i 3m T -3mA T i -3mB T i -3mC T i + 3m3m T ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ ⇒ C = 1 12E N i=1 E i • f(A i ) + f(B i ) + f(C i ) + 9 • f (A i + B i + C i ) 3 ,
where f(v) = (vm)(vm) T .

Appendix B.

In this following, a description of the theory of spherical harmonics [START_REF] Healy | FFTs for the 2-sphere-improvements and variations[END_REF] is given.

Let S 2 denote the unit sphere, which is the sphere with radius 1 and center at the coordinates' origin. Then a point v( , ) on the surface of S 2 can be expressed as v( , ) = (cos sin , sin sin , cos ), [START_REF] Hilaga | Topology matching for fully automatic similarity estimation of 3D shapes[END_REF] where 0 is measured counter-clockwise on the yz plane and 0 < 2 is measured clockwise on the xy plane. Let L 2 (S 2 ), denote the Hilbert space of square integrable complex functions in S 2 . Then the inner product of two functions f , h ∈ L 2 (S 2 ) is given by f, h =

where k l,m is a normalization constant and P m l (cos ) is the associated Legendre polynomial of degree l and order m. A visual representation of the spherical harmonics up to the 2nd degree is given in Fig. 11.

Each coefficient f (l, m) equals to the projection of the spherical function f ( , ) to the equivalent spherical

We consider band-limited functions, so the following theorem [START_REF] Healy | FFTs for the 2-sphere-improvements and variations[END_REF] to compute the SHT is used: play an analogous role with the sin factor in the integrals of Eq. [START_REF] Osada | Matching 3D models with shape distributions[END_REF].

Fig. 1 .

 1 Fig.1. The stages of the proposed 3D shape matching scheme for retrieval.

Fig. 2 .

 2 Fig. 2. (a) Initial model; (b) aligned model using the continuous principal component analysis.

Fig. 3 .

 3 Fig. 3. (a), (d) columns: alignment using NPCA; (b), (c) columns: alignment using CPCA. The objects of columns (a), (b) are better aligned using NPCA while those of columns (c), (d) are better aligned using CPCA.

  puted using CPCA and by definition, C 11 , C 22 , and C 33 represent the variance of the model's surface from the yz, xz and xy planes, respectively. By multiplying the boundaries of the spherical function with d avg in Eq. (7), the spherical functions become scaling invariant in the sense that they are scaled proportionally to the original model's size. The boundaries of the spherical functions as well as d avg and M * d avg are shown in Fig. 4.

Fig. 4 .Fig. 5 .

 45 Fig. 4. (a) Visualization of the boundaries of the spherical functions and (b) visualization of d avg and M * d avg values.

Fig. 7 .

 7 Fig. 7. Representation of a model as a set of spherical function points representing the intersections of the model's surface and occluded parts under specific viewpoints.

  [31] software setting B = 64 for the sampling bandwidth as the best choice. However this software had the limitation of computing transformations of only B = 2 n . Instead, we used the S2kit [31] software which enables transformations at any bandwidth and we set B = 40 without reducing the discriminative power. As a result, the shape descriptor's extraction time dropped significantly and the number of rays dropped from 2 * 64 * 2 * 64 = 16 384 to 2 * 40 * 2 * 40 = 6400.

Fig. 9 .

 9 Fig. 9. Precision-Recall plots for the CRSP, LF, DSR472, SH-GEDT, DB, SIL and RSH descriptors for all datasets.

Fig. 10 .

 10 Fig. 10. Example queries and the corresponding top 16 matches from the Princeton Shape Benchmark test dataset, using the concrete radialized spherical projection descriptor (CRSP) and the light field descriptor (LF). The top left model is the query model. The first row shows the results 1-8 while the second row shows the results 9-16.

Fig. 10 .

 10 Fig. 10. (continued).

(

  the sequel, we use CPCA to compute the covariance matrix C of the model's surface. If g(v) = vv T , then C is given by+ bB i + (1ab)C im) • (aA i + bB i + (1ab)C im) aA i + bB i + (1ab)C im)

Fig. 11 .

 11 Fig. 11. Representation of [Re{Y m l ( , )}] 2 for the spherical harmonics up to degree l = 2.

  harmonic function, that isf (l, m) = f ( , ), Y m l ) • e -im d • P m l (cos ) • sin d .

Theorem 1 .

 1 Let f ∈ L 2 (S 2 ) have bandwidth B. Then for each |m| l < B: j , k )e -im k P m l (cos j ), (23)where j = (2j + 1) /(4B), k = 2 k/(2B) and weights a (B) j

Table 1

 1 Cardinalities of 3D model datasets and number of classes per dataset

	3D model dataset	#Classes	#Models
	PSB test	92	907
	PSB train	90	907
	NTU	47	549
	MPEG-7	102	1300
	CCCC	55	472

Table 2 Comparison

 2

	CRSP	67.9	40.5	52.8	66.8
	LF	65.7	38.0	48.7	64.3
	SH-GEDT 55.6	30.9	41.1	58.4

of the quantitative measures for the CRSP, LF and SH-GEDT descriptors on the PSB test dataset PSB test Nearest neighbor (%) First tier (%) Second tier (%) DCG (%)

Table 3

 3 Quantitative measure scores for the proposed CRSP shape descriptor on various datasets

		Nearest neighbor (%) First tier (%) Second tier (%) DCG (%)
	PSB train 70.3		41.8		55.3	69.0
	NTU	72.1		42.0		55.3	70.1
	MPEG-7 84.8		54.9		66.2	79.6
	CCCC	81.1		54.8		68.2	78.6
	Table 4
	Quantitative measure scores for the CRSP descriptor on the "coarse2Test"
	and "coarse2Train" classifications of the Princeton Shape Benchmark in
	comparison to the DCG score of the SH-GEDT and LF descriptors
	CRSP					SH-GEDT LF
			Nearest	First Second DCG DCG	DCG
			neighbor tier	tier	(%)	(%)	(%)
			(%)	(%) (%)
	coarse2Test
	Vehicle (245)	91.8	40.5 67.0	84.9 83.1	81.8
	Animal (155)	84.5	39.0 60.5	82.3 79.6	81.2
	Household (185) 77.3	20.5 34.6	74.0 72.4	75.6
	Building (47)	34.0	13.8 22.2	54.3 50.3	59.7
	Furniture (94)	78.7	22.2 32.6	70.4 66.6	72.7
	Plant (60)	88.3	35.6 48.2	76.6 69.3	69.9
	coarse2Train
	Vehicle (230)	94.3	44.9 73.0	87.4 86.3	84.0
	Animal (123)	85.4	40.3 59.7	82.2 79.2	80.1
	Household (219) 82.6	23.9 42.9	77.4 75.5	78.5
	Building (53)	35.8	13.6 21.9	54.7 54.6	57.5
	Furniture (104)	87.5	23.3 36.1	73.0 66.8	73.5
	Plant (78)	88.5	43.7 60.6	82.3 69.4	73.6

Table 5

 5

	Comparison between the CRSP, LF and SH-GEDT descriptors in the "plant"
	class of the "coarse2Test" classification of the Princeton Shape Benchmark
	and its subclasses, using the DCG score
		CRSP	SH-GEDT	LF
		DCG (%)	DCG (%)	DCG (%)
	Plant (60)	76.6	69.3	69.9
	Subclass of plant
	Bush (9)	52.4	39.2	44.7
	Flowers (4)	58.2	34.2	35.6
	Potted plant (26)	68.1	56.3	59.9
	Barren tree (11)	71.6	47.1	47.5
	Conical tree (10)	51.0	51.2	48.2
	The cardinality of each class is given in parenthesis next to the class name.

Acknowledgment

This research was supported by the Greek Secretariat of Research and Technology (PENED "3D Graphics search and retrieval" 03 E 520).

Appendix A.

In the following a detailed description is given regarding the exact computation of a 3D model's center of mass and its surface covariance matrix, using the continuous principal component analysis [START_REF] Vranic | 3D model retrieval[END_REF][START_REF] Zaharia | 3D shape-based retrieval within the MPEG-7 framework[END_REF].

T be a point on the surface of a model; then using barycentric coordinates, v is equal to