G Rosi

I Scala

V H Nguyen

S Naili

R Vayron

G Haiat

H Yahia

S Seuret

S Jaffard

Caractérisation de la Réponse Ultrasonore d'Implant Dentaire : Simulation Numérique et Analyse des Signaux

The long-term success of a dental implant is related to the properties of the bone-implant interface. It is important to follow the evolution of bone remodeling phenomena around the implant. To date, there is no satisfactory method for tracking physiological and mechanical properties of this area, and it is difficult for clinicians to qualitatively and quantitatively assess the stability of a dental implant. In this context, methods based on ultrasound wave propagation were already successfully used by our group, in the qualitative and quantitative evaluation of primary and secondary stability of dental implants. In this study we perform numerical simulations, using the finite element method, of wave propagation in a dental implant inserted into bone. To simplify the calculations, an axisymmetric geometry is considered. Given the importance of monitoring of peri-prosthetic area, particular attention is given to the boundary conditions between the implant and the bone. The numerical results are compared with those from experimental tests. These results, numerical and experimental, are then analysed with signal processing tools based on multifractal methods. Analysis of the first results shows that these methods are potentially efficient in this case because they can explore and exploit the multi-scale structure of the signal.

Introduction

The development of acoustical modeling and of the associated numerical simulation are mandatory in order to understand the interaction between an ultrasonic wave and the bone-implant system because it can help improving the performances of the device under development. Moreover, using modelling approaches is the only solution in order to discriminate the effects of the different bone parameters (such as bone structure, geometry and material properties) on the ultrasonic response of the implant, which is nearly impossible to achieve in vivo because all parameters vary in parallel.

This study highlights the potentiality of quantitative ultrasound techniques to study dental implant stability because it gives more insight on the ultrasonic propagation in the bone-implant system [START_REF] Vayron | Variation of the ultrasonic response of a dental implant embedded in tricalcium silicate-based cement under cyclic loading[END_REF][START_REF] Vayron | Assessment of in vitro dental implant primary stability using an ultrasonic method[END_REF]. Future work should focus on considering other dental implant geometries and developing more accurate modeling of osseointegration phenomena, thus leading to a more realistic description of the medical device under development. Another perspective consists in solving the inverse problem given by the determination of the implant stability based on its ultrasonic response.

Finite Element Analysis

Due to the symmetry with respect to the implant's axis of the system and boundary conditions, an axisymmetric assumption was considered. The experimental device uses a contact planar 10MHz transducer positioned at the emerging surface of the implant [START_REF] Vayron | Variation of the ultrasonic response of a dental implant embedded in tricalcium silicate-based cement under cyclic loading[END_REF].

The geometrical configuration used in the simulation is shown in Fig. 1. The bone sample is modeled as a bi-layer structure which consists of a 1 mm thick layer of cortical bone Ω c and a trabecular bone halfspace Ω t . The geometry of the titanium dental implant Ω i (length L=11.5 mm and diameter D =4 mm) represents a dental implant commercialised by Implants Diffusion International (IDI1240, IDI, Montreuil, France). A specific healing abutment (IDI, Montreuil, France) was used and was screwed on the top of the dental implant. The healing abutment is inserted into the implant to avoid tissue formation inside the implant. The implant is completely buried in the bone sample, this configuration modelling the clinical configuration. The region Ω L is the peri-prosthetic layer, and it is defined by the domain located between the implant surface and mature bone tissue, and containing either bone tissue or water depending on the configuration. All media considered in this model are assumed to have homogeneous isotropic mechanical properties, and volume forces are neglected.

The axial symmetric equations of motion are:

ρü r - ∂σ rr ∂r - 1 r ∂σ rz ∂z - σ rr -σ θθ r = 0, ( 1
)
ρü z - ∂σ zz ∂z - σ rz r = 0, ( 2
)
where ρ is the mass density; u r and u z are the radial and axial components of the displacement vector; σ rr , σ rz , σ θθ , σ zz are the components of the stress tensor denoted by σ and the double superimposed dot stands for a second partial dime derivative. The constitutive relation using Hooke's laws is given by: σ = λtr(ε)I + µε

where λ, µ denote the Lamé constants, tr() the trace of a tensor, I is the identity tensor and ε is the strain tensor whose nonzero components are given by:

ε rr = ∂u r ∂r , ε θθ = u r r , ε rz = 1 2 ∂u r ∂z + ∂u z ∂z , (4)
ε zz = ∂u z ∂z . (5)
The emitted signal from the planar ultrasonic contact transducer is modelled as a time pulse uniform pressure applied on the top surface of the implant. Here, a timedependent function defined by

p(t) = Ae -4( f c t-1) 2 sin(2π f c t), ( 6
)
where A is an amplitude, f c is the central frequency of the pulse and t denotes the time. The discretization of the equations leads to a transient linear dynamic problem. The full description of the finite element problem can be found in [START_REF] Vayron | Finite element simulation of ultrasonic wave propagation in a dental implant for biomechanical stability assessment[END_REF]. The output radio frequency (rf ) signal was determined by computing the spatial average of the pressure at the upper surface of the implant. Extracting informations form a complex signal obtained from the ultrasonic probe is a challenging task. Indeed, the first approach used in [START_REF] Vayron | Variation of the ultrasonic response of a dental implant embedded in tricalcium silicate-based cement under cyclic loading[END_REF][START_REF] Vayron | Assessment of in vitro dental implant primary stability using an ultrasonic method[END_REF] is based on the analysis of the envelope of the signal. To verify if the numerical simulations are in agreement with experiments, the same signal processing technique as the one employed in [START_REF] Vayron | Variation of the ultrasonic response of a dental implant embedded in tricalcium silicate-based cement under cyclic loading[END_REF] was used to derive a quantitative indicator I based on the envelope of the transient signal:

I = N i=1 S (it 0 ) ( 7
)
where N is the number of samples, t 0 the sampling rate, and S (t) the envelope of the signal.

As an example, the synthetic signal recorded in the numerical tests is plotted in Fig. 3. The computation of the aforementioned indicator I is plotted in Fig. 4. As it can be seen, the numerical results are in agreement with experiments. However, since the measured signal is the result of multiple reflections occurring also at the interface between the bone and the implant, we do expect that more informations could be extracted from a multi-scale analysis.

Signal processing and multifractal analysis

With the aim of extracting more information from the complex output signal, we applied multifractal techniques to both experimental and numerical data.

Multifractal analysis aims at describing both the scaling properties of a signal and the distribution of its local singularities. The notion of "multifractal signal" was used with success in turbulence to find a formula relating the "quantity" of singularities of a given order of the signal with some structure function (which is a global quantity computed on the signal). The main interest of this formula is that, although the singularities of a signal are very complicated to localize and to estimate, the structure function is directly computable with various methods (based on increments, wavelets or Fourier transform), and the algorithms involved are relatively fast for signals of reasonable size. So the multifractal approach allows one to access to the intern structure of signals. In this framework, multifractal tools have been used to analyse the various scaling invariance properties of both the synthetic and real signals under consideration. One can hope to exhibit relevant parameters characterising the composition and/or the microstructural properties of the bone-implant (for instance, self-similar properties).

Wavelet bases

Orthonormal wavelet bases on the real set R are of the following form: There exists two functions ϕ(x) and ψ(x) such that the functions ϕ(xk), k ∈ Z, and 2 j/2 ψ(2 j xk), k ∈ Z, where Z is the set of integers, j ≥ 0, form an orthonormal wavelet basis of L 2 (R), the space of squared integrable functions. This basis is called "r-smooth" if ϕ and ψ are in the space of r-times differentiable continuous functions C r and if the derivatives ϕ (α) , and ψ (α) , for |α| ≤ r, have fast decay. This requirement implies that, for any polynomial P of degree less than r,

R P(x)ψ(x)dx = 0. ( 8
)
The orthonormal basis condition means that ∀ f ∈ L 2 (R),

f (x) = k∈Z c k ϕ(x -k) + ∞ j=0 k∈Z c j,k ψ(2 j x -k), (9)
where the c j,k and c k are called the wavelet coefficients of f and given by

c j,k = 2 j R f (x)ψ(2 j x -k)dx, (10)
c k = R f (x)ϕ(x -k)dx. ( 11
)
Note that the computation of the coefficients given by ( 10) and (11) makes sense with very little assumption on f (a wide mathematical setting is supplied by tempered distributions, if the wavelets belong to the Schwartz class: The integrals have then to be understood as duality products in the sense of distributions). We will use more compact notations for indexing wavelets:

• Instead of using the indices ( j, k), we will use dyadic intervals

λ (= λ( j, k)) = k 2 j , k + 1 2 j ,
and, accordingly:

c λ = c j,k and ψ λ (x) = ψ(2 j x -k).
Indexing by dyadic intervals will be useful in the sequel because the interval λ indicates the localization of the corresponding wavelet. Note that the wavelet ψ λ is essentially localized near λ.

• Λ j will denote the set of dyadic interval λ which index a wavelet of scale j, i.e., wavelets of the form ψ λ (x) = ψ(2 j xk).

The wavelet structure function

The wavelet structure function of f is

S f (p, j) = 2 -j λ∈Λ j |c λ | p . ( 12
)
Its order of magnitude in the limit of small scales can be related with the Sobolev regularity of the data as follows.

Let us define the wavelet scaling function by

∀p ≥ 1, η f (p) = lim inf j→+∞ log S f (p, j) log(2 -j ) . (13)
The definition supplied by (13) yields the scaling function through linear regressions in log-log plots. Let p ≥ 1 and s ∈ R; a tempered distribution f belongs to L p,s (R d ) if its Fourier transform f is a function satisfying the following property: If ĝ(ξ) = (1 + |ξ| 2 ) s/2 f (ξ), then g ∈ L p . The different variants in the definition of Sobolev spaces and the embeddings between them imply that

η f (p) = p sup s : f ∈ L p,s loc . (14)
Note that this result only holds if the wavelets used are smooth enough. The rule of thumb is that wavelets should be smoother and have more vanishing moments than the regularity index appearing in the definition of the function space. In the following we will never specify the required smoothness, and always assume that smooth enough wavelets are used (the minimal regularity required being always easy to infer).

Numerical results

A first multifractal analysis, based on the wavelet coefficients [START_REF] Abry | Irregularities and Scaling in Signal and Image Processing: Multifractal Analysis, Benoit Mandelbrot: A Life in Many Dimensions[END_REF][START_REF] Abry | A bridge between geometric measure theory and signal processing: Multifractal analysis, Operator-Related Function Theory and Time-Frequency Analysis[END_REF], shows that signals obtained from simulations and experiments have a specific particular multifractal signature. This means the log-log regressions on the multi-scale exponents could be used for a classification of the different configurations of interest. The relation (13) is used when the log-log plot regression of the wavelet structure function as the function of the scale yields a straight line. In Fig. 5 we plotted two examples of log S f (p, j) as a function ofj. The curves have been computed using the Wavelet Leader and Bootstrap based MultiFractal analysis (WLBMF) toolbox, [START_REF] Wendt | Bootstrap for Empirical Multifractal Analysis[END_REF][START_REF] Wendt | Wavelet leaders and bootstrap for multifractal analysis of images[END_REF]. First of all, this figure shows that there is no scale invariance, so that only keeping an average slope would be too reductive. On the opposite, we believe that introducing classification tools on the whole graph would be pertinent. Indeed, from the comparison between 5(a) and 5(b), we can observe that the multifractal signature for two experimental configurations enough different is relatively identical. This permanent feature points out that we must get the information from the shape of these curves. As a consequence, their classification, for a large amount of numerical and experimental data, will be the object of a further study. This study will be focused on the process of inspecting data, whose goal is discovering useful information and suggesting conclusions by classifying the signals from the shape of the wavelet structure function, as done in [START_REF] Abry | Multiscale Anisotropic Texture Analysis and Classification of Photographic Prints: Art scholarship meets image processing algorithms[END_REF].

Figure 1 :

 1 Figure 1: Cross-section view of the 3-D axisymmetric geometrical configuration used in the numerical simulations (image extracted from [4]).

Figure 2 :

 2 Figure 2: Cartography of amplitude of V z (geometrical configuration of the Fig. 1) at six different times after the emission [(a) 0.2, (b) 0.6, (c) 1, (d) 1.6, (e) 2.8, and (f) 5 µs] (image extracted from [4]).

Figure 3 :

 3 Figure 3: Numerical signal obtained from the sensor (image extracted from [4]).

Figure 4 :

 4 Figure 4: Variation of the simulated (dashed lines) and experimental (solid lines) indicator I as a function of the number of unscrewing rotation of the dental implant in bone tissue (image extracted from [4]).

Figure 5 :

 5 Figure 5: Example of log-log plots of the wavelet structure function for p = 1 : log 2 S f (p, j) versusj . The analysis is performed on experimental data obtained from two different specimens for two different stability conditions: (a) implant fully inserted in the bone, (b) implant half unscrewed and immersed into water.

Acknowledgements

This work was supported by French National Research Agency (ANR) through the PRTS program (project OsseoWave ANR-13-PRTS-0015 and by the Université Paris-Est through the PEPS program (15R03051A-METCARMAT).


Zahraa Yasseen
email: zyasseen@gmail.com

A Verroust-Blondet

A Nasri

Shape Matching By Part Alignment Using Extended Chordal Axis Transform

Keywords: Shape descriptors, shape matching, salient features, ordered part correspondence, Chordal Axis Transform, dynamic 2D shape segmentation, Dynamic Programming

One of the main challenges in shape matching is overcoming intra-class variation where objects that are conceptually similar have significant geometric dissimilarity. The key to a solution around this problem is incorporating the structure of the object in the shape descriptor which can be described by a connectivity graph customarily extracted from its skeleton. In a slightly different perspective, the structure may also be viewed as the arrangement of protruding parts along its boundary. This arrangement does not only convey the protruding part's ordering along the anti clockwise direction, but also relates these parts on different levels of detail. In this paper, we propose a shape matching method that estimates the distance between two objects by conducting a part-to-part matching analysis between their visual protruding parts. We start by a skeleton-based segmentation of the shape inspired by the Chordal Axis Transform. Then, we extract the segments that represent the protruding parts in its silhouette on varied levels of detail. Each one of these parts is described by a feature vector. A shape is thus described by the feature vectors of its parts in addition to their angular and linear proximities to each other. Using dynamic programming, our algorithm finds a minimal cost correspondence between parts. Our experimental evaluations validate the proposition that part correspondence allows conceptual matching of precisely dissimilar shapes.

Introduction

Shape descriptors and similarity measures have been a central research issue for decades. One of the main motivations is to obtain a retrieval mechanism over a database of images. Proposed approaches fall into numerous categories [START_REF] Yang | A survey of shape feature extraction techniques[END_REF] such as contour-based, region-based, transforms, moments, scale space, or spatial interrelation. Many challenges and requirements withstand shape descriptors. One with vital importance and considerable difficulty is tolerance for intra-class variety. Objects may be labeled the same whilst having disparate geometries (see Fig. 1). While the tails of the two airplanes in Fig. 1 (or the backs of the two chairs) have dissimilar shapes, they match within their contexts. Their geometric properties such as size and angle of protrusion become relative to the other protruding parts of the shape. This part alignment process allows matching conceptually similar but precisely dissimilar shapes. Although the part correspondence shown in Fig. 1 is trivial for a human, it is a composite process of dynamic segmentation, numeric representation, and search for an optimal solution when it comes to machines. The part that is concerned with dynamic segmentation also belongs to human intuition. In fact, many studies in cognitive science support the theory that humans visually segment objects into simpler parts [START_REF] Elias | Geometric determinants of shape segmentation: Tests using segment identification[END_REF][START_REF] Neri | Wholes and subparts in visual processing of human agency[END_REF]. Propositions concerning the visual segmentation process draw relations with geometric primitives [START_REF] Biederman | Recognition-by-components: A theory of human image understanding[END_REF], negative minima of curvature along the boundary [START_REF] Bertamini | Processing convexity and concavity along a 2-d contour: figureground, structural shape, and attention[END_REF][START_REF] De | Segmentation of object outlines into parts: A large-scale integrative study[END_REF], or skeleton-based segmentation [START_REF] Singh | Oxford Handbook of Perceptual Organization[END_REF] of the shape 1 . Moreover, these studies grant theories concerning numeric representations that contribute in the segmentation process. Cohen and Singh [START_REF] Elias | Geometric determinants of shape segmentation: Tests using segment identification[END_REF] conducted a series of psychophysical experiments to quantify the influences of segment length, turning angles at part boundaries, and part protrusion on segment identification. De Winter and Wagemans [START_REF] De | Segmentation of object outlines into parts: A large-scale integrative study[END_REF] analyzed the role of curvature singularities, symmetry, elongation, proximity, good continuation, and cognitive influences on rules of part formation. Bertamini and Wagemans [START_REF] Bertamini | Processing convexity and concavity along a 2-d contour: figureground, structural shape, and attention[END_REF] presented a review on evidence from behavioral, neurophysiological, imaging, and developmental studies to conclude that "convexity and concavity information along a contour is something that the visual system computes and represents as part of how shape is analyzed". Human perception is known to easily overcome imprecision and occlusion [START_REF] Rai | Fragmentary shape recognition: A bci study[END_REF]. Hence, we presume a role for these theories while attempting to match the silhouettes of objects taken from different perspectives, with part occlusion, or intra-class diversity. This role helps analyze shapes from the same perspective in which they are perceived by a human inspector: part-based representation independent of level of detail. In this paper, we propose a shape matching scheme inspired by the aforementioned theories and validate its correctness on benchmark datasets. We compute both a skeletal description and a dynamic 2D shape decomposition.

The input data is a binary image representing the silhouette of a single object. We extract the contour, locate corner points, and sample the in-between contour fragments uniformly. The advantage of locating corner points is the inclusion of the sharp features in the sample set. The region is then triangulated using constrained Delaunay triangulation (CDT). The rectified chordal axis transform [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF] (CAT) which is a descendant of the medial axis transform, and a set of pruning and merging operations provide a skeleton with an association between skeletal segments and subregions. Using this skeleton, we segment the shape and embed the segmentation in a hierarchy where the leaves are the protruding parts corresponding to the terminal segments defined by the CAT.

The subtrees in this hierarchy group terminal segments with their adjacent internal ones. After pruning away larger subtrees using some heuristics, the remaining ones represent the features of the shape on different levels of detail. We refer to all such subtrees and leaves as the visual parts that may 1 Hoffman and Richards [START_REF] Hoffman | Parts of recognition[END_REF] introduced the 'minima rule' which states that human vision defines part boundaries at negative minima of curvature on silhouettes. Biederman [START_REF] Biederman | Recognition-by-components: A theory of human image understanding[END_REF] proposed a recognition-by-components theory that assumes the existence of a finite set (total of 36) of simple geometric components that make up an object-part vocabulary. Cave and Kosslyn [START_REF] Backer Cave | The role of parts and spatial relations in object identification[END_REF] questioned whether the segmentation is performed prior to recognition and aggravate the role of spatial relations among components. Neri [START_REF] Neri | Wholes and subparts in visual processing of human agency[END_REF] presented experimental evidence to support the recognition-by-part hypothesis as opposed to recognition-as-a-whole. Delorme et al. [START_REF] Delorme | Key visual features for rapid categorization of animals in natural scenes[END_REF] studied the influence of the presence of diagnostic animal parts and their relative spatial configuration on visual processing. Higgins and Stringer [START_REF] Higgins | The role of independent motion in object segmentation in the ventral visual stream: Learning to recognise the separate parts of the body[END_REF] proposed a role for independent motion of body parts in its segmentation. contribute to the matching process.

Visual Parts are always arranged in their anti-clockwise order of appearance along the boundary of the object. Their spatial and angular distances comprise an inter-distance matrix relating every pair of them. Each visual part is described by a set of geometric attributes including relative size, eccentricity, circularity, and convexity in addition to a shape signature. The geometric attributes are multiplied by weights and assembled into a feature vector. The distance between two visual parts is then defined by the Euclidian distance between feature vectors plus the squared distance between signatures. The matching process is performed by topologically rotating the position of the 'first part' in the shape. Consequently, more acceptable segmentations of an object will result from matching it to those in the same class. For every configuration, we employ dynamic time warping (DTW) to find a minimal cost match.

Basically, the components of our work are a skeleton-based segmentation, geometric salience measures to describe visual parts and their spatial relations, and a matching scheme built on a part-to-part correspondence between shapes. The shape descriptor is invariant under affine transformations such as scaling, translation, and rotation. It has also demonstrated tolerance for occlusion, articulation, and minor deformations. The contributions of our work are as follows:

• We bring the CAT into shape description and strengthen some of its weaknesses that are revealed by dissimilarities between skeletons of objects of the same class. It is worth mentioning that Prasad [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF] does not provide a matching technique adapted to the CAT representation method.

• We alter the customary practice followed when using DTW in closed shapes by rotating visual parts instead of sample points to find the best match.

The matching process based on the protruding parts is simple and visually intuitive compared with the ones based on skeleton matching. Moreover, this matching process can compare two shapes at different levels of detail or with different topology in a unified way. The rest of the paper proceeds as follows. We discuss related work that describe shapes by parts or use their skeletons in the process in Section 2. We present our contribution to the CAT segmentation in section 3, whereas the topological structure that embeds different levels of details is discussed in Section 4. The geometric properties, saliency, and distance measure of parts are detailed in section 5. In section 6, we describe the DTW method, and present the retrieval results obtained on Kimia-99 and Kimia-216 datasets [START_REF] Sebastian | Recognition of shapes by editing their shock graphs[END_REF]. Section 7 presents a brief conclusion of the work presented in this paper.

Related Work

Hoffman and Richards [START_REF] Hoffman | Parts of recognition[END_REF] introduced the 'minima rule' which states that human vision defines part boundaries at negative minima of curvature on silhouettes. Following the 'minima rule' [START_REF] Hoffman | Parts of recognition[END_REF], most of the fragment-based shape descriptors [START_REF] Mohammad | Classification of silhouettes using contour fragments[END_REF][START_REF] Longin | Shape similarity and visual parts[END_REF][START_REF] Ma | From partial shape matching through local deformation to robust global shape similarity for object detection[END_REF][START_REF] Kang | Classification of contour shapes using class segment sets[END_REF][START_REF] Wang | Bag of contour fragments for robust shape classification[END_REF] extract boundary segments by locating critical points based on convexity or curvature variations. Latecki et al. [START_REF] Longin | Shape similarity and visual parts[END_REF] partition the contour of a shape into an ordered set of convex and concave boundary arcs. In a somewhat similar way to our approach, they assume that a single visual part of one curve can correspond to a sequence of consecutive parts of the second curve and employ dynamic programming to find a minimal cost correspondence. Sun and Super [START_REF] Kang | Classification of contour shapes using class segment sets[END_REF] extract redundant and overlapping contour segments between all the pairs of critical points at extremal curvature values. They represent the segments in a vector space where, after applying PCA, they note that the first few basis 'approximately capture height, sideways shearing, diagonal stretching, and inflation/deflation of a segment.' Daliri et al. [START_REF] Mohammad | Classification of silhouettes using contour fragments[END_REF] extract critical points from the local extrema of the scale invariant curvature profile and construct a dictionary of 11 symbols representing types of contour fragments. The symbols signify curve parts, straight lines, sharp angles, concave and convex parts and such. Guo et al. [START_REF] Guo | Computing importance of 2d contour parts by reconstructability[END_REF] introduced the notion of 'part importance' based on reconstruction quality, variation, and uniqueness and measured by some conditional entropy. Ma and Latecki [START_REF] Ma | From partial shape matching through local deformation to robust global shape similarity for object detection[END_REF] propose a shape descriptor with basic geometric units similar to the shape context [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF] and suitable for matching of edge fragments. Applying discrete contour evolution [START_REF] Bai | Skeleton pruning by contour partitioning with discrete curve evolution[END_REF], Wang et al. [START_REF] Wang | Bag of contour fragments for robust shape classification[END_REF] locate critical points to fragment the shape's contour. They use the bag of words model to represent a shape described by 400 overlapping fragments with a shape context descriptor of 300 bins each. Laiche et al. [START_REF] Laiche | Curve normalization for shape retrieval[END_REF] explore the least squares model to represent a shape by a set of normalized cubic curves.

Skeleton-based shape descriptors either require a tree or graph matching method [START_REF] Erdem | A similarity-based approach for shape classification using aslan skeletons[END_REF][START_REF] Sebastian | Recognition of shapes by editing shock graphs[END_REF][START_REF] Sebastian | Shock-based indexing into large shape databases[END_REF][START_REF] Siddiqi | Shock graphs and shape matching[END_REF][START_REF] Sundar | Skeleton based shape matching and retrieval[END_REF][START_REF] Torsello | A skeletal measure of 2d shape similarity[END_REF][START_REF] Torsello | Discovering shape classes using tree edit-distance and pairwise clustering[END_REF] or escape the incurred instabilities by just drawing on topological characteristics provided by the skeleton. Instead of directly matching skeletons, Xie et al. [START_REF] Xie | Shape matching and modeling using skeletal context[END_REF] use the MAT to find an optimal point correspondence between two shapes. Goh [START_REF] Goh | Strategies for shape matching using skeletons[END_REF] presents a discussion on strategies that should be adopted in skeleton-based shape matching. Namely, he stresses on the importance of dynamic part-decomposition, visual saliency measures, and the use of both global and local measures. Aslan et al. [START_REF] Cagri Aslan | Disconnected skeleton: Shape at its absolute scale[END_REF] represented symmetry branches by vectors from an origin located by an evolution process to the disconnection points where the limbs connect to the body. Supported by perceptual motivation, Temlyakov et al. [START_REF] Temlyakov | Two perceptually motivated strategies for shape classification[END_REF] propose strategies to improve shape classification. They identify thin and elongated strand structures and their base structure and measure similarity between shapes on these structures separately. Feldman et al. [START_REF] Feldman | An Integrated Bayesian Approach to Shape Representation and Perceptual Organization[END_REF] propose a skeleton that implies part decomposition. Similarity is estimated by the probability that one shape would 'grow' from the skeleton of the other. Lin et al. [START_REF] Lin | Shape classification using hybrid regional and global descriptor[END_REF] propose a hybrid shape descriptor that draws on skeletal and contour properties of the shape.

Bai et al. [START_REF] Bai | Skeleton-based shape classification using path similarity[END_REF] propose a matching method that views the skeleton as a set of shortest paths between pairs of end points. They apply a skeleton pruning process on the medial axis based on discrete curve evolution [START_REF] Bai | Skeleton pruning by contour partitioning with discrete curve evolution[END_REF]. In a subsequent publication [START_REF] Bai | Integrating contour and skeleton for shape classification[END_REF], the skeleton paths shape descriptor is used in conjunction with contour fragments to design a shape classifier. Using a similar endpoints representation for the shape skeleton, Giang et al. [START_REF] Ngo Truong | Skeleton-based shape matching using higher-order constraints[END_REF] further compare tuples of terminal nodes. They define a hyper-graph that relates a node to its k-nearest neighbors.

It is quite evident that the skeleton embeds crucial information about the shape. The above methods provide a rich set of ideas in this regard. In addition to Goh's elaborate study [START_REF] Goh | Strategies for shape matching using skeletons[END_REF] on important features that should be incorporated in any matching method, it seems that there is a general consensus on the distinguished role of endpoints of the skeleton or the terminal parts of the shape.

Skeletonisation and Segmentation

The skeleton that we seek is different from the medial axis transform (MAT) of a shape. This skeleton is not only used to depict the axial representation of the shape, but must also provide an intuitive segmentation. In other words, we seek a skeleton-based segmentation where the main visually protruding parts correspond to regions or to union of regions of this segmentation. Although the MAT does not provide such a part-to-axial segment correspondence, the CAT, which is based on principles borrowed from the MAT construction, provides a skeleton and an associated segmentation conjointly. In this section, we give a brief explanation on how the CAT segments and skeletonizes 2D shapes. We also shed some light on inadequacies of the CAT in its current form [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF] and propose various amendments to qualify it for shape description and matching.

The CAT in Its Present Form

The medial axis transform (MAT) of a 2D shape bounded by its contour line C is the locus of the centers of disks interior to the shape and tangent to C in two or more points. The formation of the chordal axis transform (CAT) is based on the observation that the maximal discs used in MAT identify the triangles of the constrained Delaunay triangulation (CDT) of the shape's interior. In a CDT, all vertices lie on the boundary of the shape, giving rise to three types of triangles: terminal, sleeve, and junction (see Fig 2). A terminal triangle has two of its edges on the shape's boundary in addition to an interior edge. Sleeve and junction triangles have two and three internal edges respectively. The term chord is used to refer to internal or non-boundary edges in the CDT.

In its primitive form, the CAT is constructed by joining the midpoints of the chords and the centroids of junction and terminal triangles of the CDT. In addition to the lack of an evident segmentation, the main problem with this version of the CAT is the exclusive third degree topology of junction points located at the centroids of junction triangles. The rectified CAT [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF] approaches a solution to this problem by using polygonal elements instead of triangular elements. These polygons are defined by a subset of the chords of the CDT, denoted by the chords strength profile (CSP), and are classified into three types (see Fig 2). Terminal and sleeve polygons have one and two chords respectively amongst their edges. Junction polygons have three or more chords where the number of chords is equal to the degree of the junction point located at the centroid of the polygon. The skeleton is formed by the segments connecting mid points of the chords in addition to the centroids of junction and terminal polygons. In addition to the skeleton, these polygonal elements provide a segmentation of the shape. In what follows we describe how the subset of chords, the CSP, is selected in order to construct the skeleton and a segmentation of the shape. From here on, we will use the term "segment" or "region" to refer to the polygonal elements.

Prasad's [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF] idea for selecting the CSP is based on quantifying the notion of the strength of chords. A chord's strength is measured by estimating how far the quadrilateral formed by its flanking triangles is from being cyclic. The chords in the triangulation are ordered in an anti-clock wise traversal of the boundary. A local maximum in this set is a chord whose strength is greater than both of its neighbors and strictly greater than one of them. These local maxima make the CSP. Once the CSP is formed, the segmentation and the skeleton follow trivially (see Figs. 3(a), 4(a), and 5(a)).

Extended Chordal Axis Transform

Employing the CAT in shape representation and matching incurs additional requirements of its segmentation and skeletonisation. It also reveals some weaknesses exposed by dissimilar topologies of similar objects. Trivially, it is necessary to have tolerance for an extent of boundary variation where similar objects tend to have similar skeletons. Moreover, some objects are visually segmented into a centralized structure with one junction segment connected to most of the terminal parts (see Fig. 3(a)). Due to its reliance on the CDT, the CAT tends to preclude such junction segments favoring a distributed topology. As we want to keep only significant protruding parts, excessive branching on terminal segments should be avoided and visually unimportant segments should be merged. Nonetheless, in certain instances, terminal segments recede and adjacent junction segments take parts of what visually belong to protruding parts (see tail and fin in Fig. 4(a)).

In general, the qualifying properties for a shape descriptor are the invariance to scale, translation, and rotation and tolerance for occlusion, articulation and minor deformations. The CAT is evidently invariant to the affine transformations. However, the second batch of properties is more difficult to accomplish. In our version of the CAT, we tweak its formation phases so as to reduce its sensitivity to minor deformations. Figs. 3, 4, and5 show multiple examples of objects from the same classes (star, fish, and hand ) illustrating the difference between the CAT in [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF] and the CAT that we construct in (a) The CSP (the chords of the CDT that are highlighted in red) as proposed in [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF] (first row) and the consequent skeleton and segmentation (second row).

(b) The expanded CSP where added chords are the compactness chords shown in green and the enclosing disk chords in blue (first row) and the consequent skeleton and segmentation (second row). (a) The CSP (the chords of the CDT that are highlighted in red) as proposed in [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF] (first row) and the consequent skeleton and segmentation (second row). The CAT after expanding the CSP and applying the pruning, merging, and trimming operations results in more restrained junction areas than that proposed in [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF] (a) The CSP (the chords of the CDT that are highlighted in red) as proposed in [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF] (first row) and the consequent skeleton and segmentation (second row).

(b) The expanded CSP where added chords are the compactness chords shown in green and the enclosing disk chords in blue (first row) and the consequent skeleton and segmentation (second row).

(c) The skeleton and segmentation after applying the pruning, merging, and trimming operations on the segments shown in (b).

Figure 5: The CAT proposed in this paper results in closer skeletal structures for similar objects than the CAT in [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF]. Despite the sever part occlusion in the left most object, its CAT retained most of the structure of the un-occluded object.

Figure 6: From left to right: chords selected into the CSP in red, the concave region in dark red and the salient parts with round geometry in dark green, the chords with maximal measure of compactness in green and maximal enclosing disk in blue, and the resultant segmentation.

this paper. The figures also show the reduced difference between instances of the modified version of the CAT in the same class of shapes. In particular, the structure of the CAT withstands the sever part occlusion in Fig. 5(c). In the original settings (see Fig. 5(a)) every instance of the CAT relates the four fingers of the hand in a different topology. The previously mentioned shortcomings of CAT can be addressed by interjecting the CAT formation process in two main phases:

1. In the first phase, we expand the CSP to include compactness chords and the enclosing disk chords that will be introduced shortly (see Figs. The two phases added to the CAT formation process allow creating an intermediate finer segmentation that undergoes conditional merging operations. In what follows, we will describe these procedures in detail.

Extended CSP

In certain cases, some visually important chords are not selected into the CSP. This may result in undesirably concave segments (see the dark red segments of the bug in Fig. 6) or neglected salient parts with round geometry (see the dark green segments of the instrument in Fig. 6). The disadvantage of such cases is amplified by the vital role of the visual protruding parts in the shape's description. It is important to associate all visually protruding parts to a segment or a union of adjacent segments. The concave junction segment in the bug and the bulging sleeve of the instrument (Fig. 6) occupy parts of salient visual features. Such cases motivate searching for criteria to segment the shape other than the strength of the chord introduced in [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF].

The strength of a chord is an estimation of how the associated quadrilateral is cyclic. We incorporate the inverse of its circularity ratio denoted as measure of compactness in the CSP selection. In effect, we added more search criteria and expand the CSP to include local maxima of the following values:

Maximal compactness:

The measure of compactness of the quadrilateral formed by its flanking triangles as defined in [START_REF] Ballard | Computer Vision[END_REF]. This is the square of its perimeter divided by its area.

Maximal enclosing disk:

The chord splits the boundary of the shape into two parts. Let the chord's boundary segment be the smaller part. The value that we seek here is the number of points on the chord's boundary segment such that all of them are interior to the disk formed by rotating the chord around its midpoint. Having one point outside this disk sets this value to zero (see Fig. 7).

Having established the CSP, the shape is easily segmented by the chords belonging to this set (see Figs.

Insignificant Terminals

Generally speaking, terminal segments represent the parts protruding from the shape with a certain degree of saliency. If a terminal segment, adjacent to a junction segment, does not have sufficiently sharp diversions represented by negatively oriented angles around its boundary, it is tested for insignificance (see Fig. 8). Such a terminal is considered insignificant if the ratio of the length of its chord to the length of its boundary segment exceeds values relative to its size as follows:

1. The area of the terminal segment is smaller than ϵ and the chord length ratio is greater than µ 0 . 2. The chord length ratio is greater than µ 1 .

where values of ϵ, µ 0 , and µ 1 are 1%, 0.75, and 0.9 respectively. Terminal segments with such properties have a flat-like geometry and are thus merged to their neighboring junction segments.

Junction Blocks

Junction segments determine the topology of the skeleton. Whether this topology is centralized or contains symmetrical parts depends on the degrees and connectivity of its junctions. In general, shapes with rotational symmetry have a central junction segment connecting the rather similar visual protruding parts. Nevertheless, the skeleton of a shape with reflectional symmetry is expected to carry the same symmetry as well where mirrored parts have symmetric skeletons. These requirements are not always met by CAT in its original form [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF]. For example, the skeleton of the leaf in Fig. 9(a) is deviated from its axis of symmetry by the quasi central sixth degree junction. As such, this junction should be split in a way that allows more symmetry in the skeleton (see Fig. 9(c)). On the other hand, the octopus in Fig. 9(a) has a partial rotational symmetry where approximately similar protruding parts connect to a central unit. Thus, the block of adjacent neighboring junction segments is better merged to one higher degree junction producing a skeleton with a centralized topology (see Fig. 9(c)). The CAT presented in this paper handles these cases by expanding the CSP thus refining the segmentation and then traversing blocks of adjacent junction segments (see Fig. 9(b)) to test possible merge operations. In this section we explain the merging process in detail.

The relation between a chord

-→ C and a junction point J can be described by the length of the skeletal segment ∥ -→ S ∥ connecting J and the mid point of -→ C in addition to θ, the angle between -→ C and -→ S . A small θ or a long -→ S are indications of weakly related chord and junction point. In a block of adjacent junctions, the properties that need to be balanced out for all chords and junction points are the following:

1. Minimum difference between θ and π/2.

Minimal skeletal segment length ∥

-→ S ∥. This property prioritizes relating chords to proximal junction points. 3. Minimal junction segments degree.

Minimal chord lengths

∥ - → C ∥.
In a block of adjacent junctions, an optimal setting for the first two properties can always be reached by merging all junctions to one. For example, the largest block of junctions in the leaf in Fig. 9(b) would be merged to one degree 10 junction segment since the collected values of ten skeletal segments and ten angles will definitely be less than that of 24 skeletal segments an 24 angles. To prevent forming this large degree junction segment, we infuse the third condition that is to minimize the degree of produced junction segments. On the other hand, if the block contains a relatively large chord such as the highlighted chords in the octopus of Fig. 9(b), it means that merging its two junctions will produce a junction segment that has less disparate dimensions. For this reason and to balance out the tendency towards smaller junction degrees created by the third condition, we also incorporate the chords' lengths in the fourth condition. To formulate a numerical representation for the above properties, we define the strain of a junction segment. Let n be the number of chords of a junction segments and θ i the angle between the i th chord -→ C i and the i th skeletal segment -→ S i . The strain of this junction is quantified by the following formula:

n × n ∑ i=1 (θ i - π 2 ) 2 × ∥ -→ C i ∥ × ∥ - → S i ∥
Since these blocks are normally sets with small cardinalities, it is easy to test all possible combinations of merging. The set that we consider for every attached block of junctions is that of their (a) The CAT as proposed in [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF] may fail to depict important shape properties such symmetry and centr.

(b) The segmentation resulting form the extended CSP produces junction segments with smaller degrees. The junction blocks that are considered for different configurations by potential merge operations are 3 blocks for the leaf (with 7, 3, and two junction segments in each block respectively) and one block in the octopus (with 3 junction segments).

(c) The optimal setting for the leaf sustains all smaller degree junctions whereas in the octopus it is achieved by merging the whole block into one 8 degree junction segment. connecting chords (see highlighted chords in Fig. 9(b)). Each chord has the possibility to be included in the final configuration. If it is not included, it means that the segments on either of its sides are merged to one. We generate all possible configurations, compute the total strain measure for each, and select the one with minimal strain. While the minimal strain setting in the leaf is achieved by sustaining all chords, it requires merging the junction block in the octopus to one 8 degree junction segment (see Fig. 9(b)).

Furthermore, since the role of a junction is more topological than geometric, we try to reduce its size by trimming away triangles from the CDT that, if removed, neither change the topology nor produce more strain.

The Topological Structure of CAT Segments

Although insignificant terminal parts are pruned, the remaining ones still differ in their level of importance as to what identifies a shape. In Fig. 10 we show similar objects that have different sets of terminal segments. A part-based matching in such cases is expected to match the wing terminal segment in one object to the set of all segments that make up the wing in the other. In order to allow this, without losing detailed features when present, we seek a hierarchical topology where levels correspond to levels of detail.

We first describe the topological structure and connectivity of the segments produced by the CAT skeletonisation. The graphical structure is evident (see Fig. 11). For each segment we associate a node in the graph where terminal, sleeve, and junction nodes have, one, two, and three or more neighbors respectively. Since our data is comprised of silhouettes without holes, the graph connecting these nodes cannot have any cycles.

The relation between features across different levels of detail within the same object is to a certain extent, a node/subgraph relation. For example, while the feathers on the wing in Fig. 10 are terminal nodes, the wing can be represented as a subgraph that includes these nodes. However, not any subgraph is a feature in a higher level of detail representation of the object. A tree structure of the CAT nodes may resolve this issue provided that the hierarchy is based on geometric as well as topological properties. It is quite obvious that the leaves of a tree representation are the terminal nodes. However, tracing the leaves back to a root is not quite as simple. Our scheme is based on simulated merge operations involving the less salient terminals with their neighborhoods in a recursive manner. At every merging stage, a new set of compound terminals emerge. These terminals are comprised of nodes constituting subgraphs that represent potential features.

Finding the root

The tree construction is based on the node depth computation. Since the node connectivity is already there, it suffices to find a root and the tree structure will follow. The root is expected to have a central position in the shape. So, in addition to the topology of the nodes, their geometry has to be incorporated in the depth computation process. We start from the terminal nodes and assign them zero height values. In an iterative process, smaller terminals are merged with their connecting nodes. For example, in the first iteration, all terminals with area and perimeter values less than σ, a percentage of the whole shape, are merged with their neighbors. Consequently larger terminals are formed and the next iteration merges terminals smaller than σ × 2. When all the terminal neighbors of an interior node (sleeve or junction) are merged, it becomes a terminal if it has only one non-terminal neighbor. In this case, it is assigned a height value equal to the index of the current iteration. The remaining one interior node is the root (see Algorithm 1). The bottom up tree shows the iterative process that assigns different levels of depth to junction nodes. The wings (compound terminals) build up as the process proceeds starting at the terminals towards the sought root. The top-down tree shows the hierarchical structure of the CAT nodes.

Data: sigma

This virtual merge operation yields a new entity which we will use the term 'wing' to refer to. The perimeter and area threshold is increased at every iteration and the remaining interior nodes are tested for merging again (see Fig. 12). The last remaining junction note is considered the root.

Visual parts' representation in the tree

Visual parts, i.e. the distinctive features in the silhouette of a shape, are not strictly the terminal nodes. This is a consequence of the level of detail of the shape at hand where terminals are the visual parts in the finest level. In Fig. 13 we show two similar objects but which are represented by trees with different topologies. The figure shows two differently sized and placed subtrees that abstractly represent the same part in these objects. So, to have all possible visual parts contribute in the matching process, we need to include some subtrees in addition to the leaves.

However, the size and complexity of the visual parts must be limited so as to preserve the partbased characteristic of the method. To that effect, and as a heuristic measure, we avoid having more than two salient features in a visual part. Salience is a visual property related to the part's size and eccentricity. We will introduce the set of conditions that determine the salience of a visual part in the next section after we introduce the geometric parameters describing it.

Each subtree viable for representing a visual part is collected in a new entity that we call a wing node (see Fig. 14). The wing construction process resembles that which finds the root. Subtrees with smaller sizes are tested first. If the number of salient visual parts that a subtree contains is equal to or less than 2, a new wing node is formed and added to the set of visual parts representing the shape.

The Shape Descriptor

The shape is described by its sequence of visual parts represented as terminal or wing nodes. Each visual part occupies a fragment in the boundary of the silhouette. It is described by a set of shape parameters expressing visual characteristics such as size, degree of protrusion, and eccentricity. Some of these parameters are only derived from the geometry of the terminal parts whereas others are measures taken relative to the whole shape. In addition to the set of numeric values computed for visual parts are:

1. The percentage of its area. 2. The percentage of perimeter it occupies along the boundary. 3. The eccentricity is the ratio of its height away from the shape to the ratio of its width along the shape. In order to compute this, we obtain the minimal bounding box using the rotating calipers algorithm around its convex hull. Then we project the edges of the shape segment on the two sides of the minimal rectangle. The side with the greater sum of projections is the one selected to be aligned with the shape. (see Fig. 15) 4. The circularity defined by 4π × p 2 /a, p and a being perimeter and area respectively. 5. Rectangularity is the ratio of the area of the segment to that of its minimal bounding rectangle. 6. Convexity is the ratio of the perimeter of the convex hull to that of the shape segment. The more convex the shape the higher its convexity. To give a greater weight to objects with lower convexity, we reversed the formula as follows:

Convexity (reversed) = (1 -Convexity) × N shape /N hull
where N shape is the number of vertices of the node and N hull is the number of vertices of its convex hull polygon. 7. Solidity is the ratio of the area of the shape segment to that of the convex hull. Similar to convexity: The depth of the interior nodes starting from 0 at the root and increasing as nodes get topologically away from it. Interior nodes are arranged according to depth which is color coded from zero to 6 (horizontal legend). Terminals are color coded according to their clockwise order (vertical legend). The yellow boxes mark subtrees corresponding to equivalent parts in different objects. The minimum bounding boxes of terminal nodes. Note that terminals A and B have similar bounding boxes, but the eccentricity of A is less than 1 whereas that of B is closer to 3.

Solidity (reversed) = (1 -Solidity) × N shape /N hull

8.

Bending energy is the sum of curvature values of each point on the boundary segment of this node. To compute the curvature, the whole object is scaled by the reciprocal of the maximal distance to the objects centroid. The discrete curvature at a boundary point B preceded by A and followed by C is:

|(B -A) × (C -B)| |B -A||C -B||A -C| 9.
The chord length ratio (CLR) is the ratio of the length of the chord connecting the node to the main object to the perimeter of its boundary segment. 10. The radial distance signature of re-sampled points. We re-sample the boundary parts of polygonal segments into N s points and compute the distance to the segments centroid. We then divide this signature by the unit length which is the ratio of the perimeter of the segment by its number of points.

In addition to their role in part to part dissimilarity measure, these parameters serve two other purposes. They are used to estimate visual saliency in the wing nodes construction phase. They also constitute the penalty added to the matching cost when the associated visual part is excluded from the correspondence. In what follows we detail out the formulae that evaluate the distance between two visual parts, the penalty, and visual salience assessment.

The distance between visual parts

Each terminal or wing node is described by two vectors: the assembly of the parameters p and the radial distance signature r. The distance between any two nodes is given by the following formula:

d(v 1 , v 2 ) = sqrt [ CLR ∑ I=AREA (p 1 [I] -p 2 [I]) 2 ] + Ns ∑ i=1 (r 1 [i] -r 2 [i]) 2 (1)
Similarly, the penalty of excluding a node is given by: sqrt

[ CLR ∑ I=AREA p[I] 2 ] + Ns ∑ i=1 r[i] 2 (2)
The visual parts are not independent of their context in the object. To represent this context, an inter-distance matrix (iDM ) saves proximity and angular relations between nodes. Every entry (i, j) in this matrix is a 3 dimensional vector where:

1. The first dimension is the ratio of the Euclidian distance between the end point of node i and the start point of node j to the perimeter of the object. 2. The second dimension is the percentage of the bending energy along the boundary segment between nodes i and j. 3. For every node, we compute the direction vector which is the average of the vectors connecting the midpoint of its chord to the points on the boundary. The third dimension is the angle between the direction vector of node i and node j.

The contribution of inter-distances in the computation of the cost of the match between two shapes involves 4 nodes. For every consecutive pairs of visual parts in the correspondence between shapes A and B, (v A i , v B s ) and (v A j , v B t ), the squared Euclidian distance between iDM A (i, j) and iDM B (s, t) is added to the total cost.

Part salience

During the wing node formation process, if the subtree contains more than two nodes with any salient feature property, the subtree is discarded. We assess the saliency of a visual part by the following tests: 1. Relatively large size: Area or perimeter percentages greater than a threshold th large 2. Relatively long shaped terminals: For perimeter values smaller than th large , if multiplying the perimeter by the eccentricity of the visual node is larger than k long , then it is considered long node. 3. Bent terminals: We estimate this for perimeter greater than th bent and convexity less than c bent . If the connected pockets of void areas created between the convex hull's perimeter and the terminal's perimeter have areas larger than th large (see Fig. 16), we consider that the feature is salient enough.

Matching Shapes

We employ dynamic time warping [START_REF] Myers | Performance tradeoffs in dynamic time warping algorithms for isolated word recognition[END_REF] to find a minimal distance correspondence between parts of matched shapes. Roughly, the idea is to alternatively rotate one shape and fix the other while calculating a distance matrix for every obtained alignment (Algorithm 2). The rotation is topological, i.e. it is performed by shifting the starting point of the terminal parts. The valid start points are terminals that are the first to appear in constituting wing nodes along the counter clockwise direction (Algorithm 3). The number of possible start points in an object is analogous to the number of skeletal branches connected to its root. Each row in the matrix is the distance between a terminal in the first shape and all terminals in the other. A minimal distance path is calculated for every matrix resulting in one-to-one, one-to-many, many-to-one, or many-to-many pairing. The matrix that produces the minimal distance among others represents the best alignment. We show some examples in Fig. 17 of the first 11 retrievals and the associated correspondence of terminal segments in each match.

Parameters

The 2D shape databases that we use are sets of binary images. We extract each object's contour and sample it to 100 points for triangulation. Adaptive sampling where the number of sample points is deducted from the shape as in [START_REF] Dilip K Prasad | A novel framework for making dominant point detection methods non-parametric[END_REF] can be applied. However, since the subsequent triangulation will affect skeleton generation, a uniform sampling of the boundary is more preferable. The visual parts' radial distance signature is taken on 12 uniformly sampled points on each part. The geometric attributes presented in section 4 have varying ranges of values. While the area and perimeter are percentage values , the convexity or solidity are less than 1. This gives much more weight in the distance computation to larger valued parameters. To amend this, we apply weights and balance their values into more reasonable ranges. For the same reason, we apply weights to the energy and angle components (second two dimensions) of the vectors of the iDM . These weights are shown in Table 1. Some parameters mentioned through out the paper are listed in Table 1. These values were deducted by inspecting the dataset where terminals with areas less than ϵ (1%) look insignificant and those with areas greater then th large (20%) of the shape look relatively large. As will be shown in the next section, these parameters do not affect the retrieval results as much as they affect the timing performance.

Complexity

The feature vector of each terminal or wing node is made of 22 numeric values representing the parameters and the signature. Table 3 shows statistics about the average number of terminals and wings per object using both the original CAT [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF] and the extended CAT on the datasets that we used. The object is also described by the inter-distance matrix which is a 3 dimensional matrix of size: (N T + N W ) × (N T + N W ) × 3 where N T and N W are the numbers of terminal and wing nodes respectively.

The matching algorithm constructs a matrix for every start point in the two matched objects. The dimensions of each matrix is N 1 T × N 2 T where N 1 T and N 2 T are the numbers of terminal nodes of the objects. Every entry in the matrix is responsible for computing the cost of all possible matches relating the associated terminals. The number of options at each entry is analogous to the number of wings the terminals belong to. Table 3 shows the average number of wings of the terminals belonging to all the objects in the testing datasets. Having said that, the total number of node distance computations as specified in Eq. 1 is:

(N 1 start + N 2 start ) × (N 1 T × N 2 T ) × (N i W oT × N j W oT )
where N start is the number of start point of an object and N i W oT is the number of wings associated with the i th terminal. As such, the time complexity of the whole algorithm is O

(N 2 start × N 2 T × N 2 W oT
). The number of wings of a terminal N W oT is regulated by the condition on the maximal number of salient parts allowed in a wing. Table 3 shows how the average values of N W oT in both datasets and for both original and extended CAT are roughly the same. This is an indication that we may consider the N 2 W oT ) factor in the time complexity notation as a constant. N start is analogous to the number of major parts of a shape. It is class related and independent of the level of detail of the silhouette. On the other hand, the finer the level of detail the greater N T becomes. Since all objects in our system are resampled to a fixed number of points, even N T becomes under control. Moreover, the average number of start points N start and terminals N T in the extended CAT and those of the original CAT [START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF] (see Table 3) show that our pruning and merging operations achieved a more than five times faster average retrieval time for both dataset.

Results

We tested our algorithm on two benchmark databases: Kimia-99 and Kimia-216 [START_REF] Sebastian | Recognition of shapes by editing their shock graphs[END_REF]. The percentages obtained for the datasets are 98.18 % and 95.83% for the nearest 10 and 11 retrievals respectively. To show the effectiveness of our algorithm in the presence of occlusion, we selected examples from Kimia-99 well known for such cases and showed their retrieval results in Fig. 18. Comparisons of retrieval rates for algorithms tested on these sets are shown in tables 4 and 5. Our method is positioned first on Kimia-216 and obtains the second best score on Kimia-99 dataset. The results show that our method outperforms the skeleton-based methods proposed in [START_REF] Bai | Skeleton-based shape classification using path similarity[END_REF][START_REF] Goh | Strategies for shape matching using skeletons[END_REF][START_REF] Sebastian | Recognition of shapes by editing their shock graphs[END_REF][START_REF] Temlyakov | Two perceptually motivated strategies for shape classification[END_REF]. In addition, there is considerable improvement over using the segmentation of the original CAT.

The effectiveness of each part of the feature vector in the distance estimation can be made sense of by a few experiments. We first set all parameter weights to zero. The average precision for Kimia-99 and Kimia-216 are 92.5% and 90.5% respectively. When all weights are set to 1, the obtained results are 93.7% and 90.6%. Turning the computations around by discarding the radial distance signature from the feature vector and resetting the weights to those in Table 1, the results drop to 90.7% and 86.9% for the two datasets. The weights of the geometric parameters (eccentricity, circularity, bending energy, rectangularity, convexity, solidity, and chord length ratio) are tuned to balance their effect with respect to each other and improve the retrieval results. Setting the weight of each parameter to zero in turn, we rerun the experiments to obtain the average precision for each dataset. The 9 precision values for Kimia-99 and Kimia-216 are in a close range around 95% and 92% respectively. Moreover, the precision recall plots for the 9 runs show how the parameters with their current weights have a somewhat similar effect on the retrieval results.

Conclusion and Future Work

We propose a new shape descriptor inspired by human vision and cognition. Our descriptor is comprised of the arrangement of visually salient parts along its boundary. Based on an extended version of the chordal axis transform, we segment the shape and embed the segmentation in a hierarchical structure. This allows capturing parts on different levels of detail. Based on a part representation by radial distance signature, and using a set of simple geometric properties, we devised a dynamic programming approach to find a minimal cost correspondence between two shapes. We show that a redundance-prone representation of overlapping boundary fragments can be avoided. Balanced on

Figure 1 :

 1 Figure 1: Matching shapes by aligning their visual parts allows corresponding parts that are conceptually similar but precisely dissimilar.

Figure 2 :

 2 Figure 2:Left to right: the constrained delaunay triangulation (CDT), junction (red), sleeve (green), and terminal (blue) triangles, the chords in the strength profile shown in a darker color, and the segmentation of the shape to three types of segments: junction (red), sleeve (green), and terminal (blue).

  (c) The skeleton and segmentation after applying the pruning, merging, and trimming operations on the segments shown in (b).

Figure 3 :

 3 Figure 3: The CAT of objects with centralized skeletal topology.

  (b) The expanded CSP where added chords are the compactness chords shown in green and the enclosing disk chords in blue (first row) and the consequent skeleton and segmentation (second row).

  (c) The skeleton and segmentation after applying the pruning, merging, and trimming operations on the segments shown in (b).

Figure 4 :

 4 Figure4: The CAT after expanding the CSP and applying the pruning, merging, and trimming operations results in more restrained junction areas than that proposed in[START_REF] Prasad | Rectification of the chordal axis transform skeleton and criteria for shape decomposition[END_REF]

  3(b), 4(b), and 5(b)). 2. In the second phase, we prune insignificant terminals and analyze blocks of adjacent junction segments for possible merging into fewer higher degree junction segments (see Figs. 3(c), 4(c), and 5(c)).

Figure 7 :

 7 Figure 7: The chord colored in green is a local maximum since its neighbors either have a smaller number of boundary points belonging to their half disks (blue chord) or have some points on their boundary segment outside their half disk (red chord).

  3(b), 4(b), and 5(b)).

Figure 8 :

 8 Figure 8: The result of merging insignificant terminal segments shown on two objects belonging to the same class. The topological dissimilarity between the two objects decreased significantly after the merging.

Figure 9 :

 9 Figure 9: Analyzing blocks of junctions produced by the extended CSP for possible merge operations. A comparison between 9(a) and 9(c) show how the skeleton produced by our method is more successful in reflecting shape properties such as symmetry.

Figure 10 :

 10 Figure 10: Similar shapes with different topologies and terminal segments.

Figure 11 :

 11 Figure 11: The topology of the CAT segments resulting in a graph structure.

Algorithm 1 :

 1 Result: Only one internal node: the root. while number of interior nodes > 1 do forall the interior nodes N do forall the terminal neighbors T do if areaOfT < sigma AND perimeterOfT < sigma then merge T into N; NeighborsOfN = NeighborsOfN -1; end end if NeighborsOfN = 1 then change N into a terminal; end end sigma = sigma × 2 end Locating the root of the tree.

Figure 12 :

 12 Figure 12:The bottom up tree shows the iterative process that assigns different levels of depth to junction nodes. The wings (compound terminals) build up as the process proceeds starting at the terminals towards the sought root. The top-down tree shows the hierarchical structure of the CAT nodes.

Figure 13 :

 13 Figure 13: The depth of the interior nodes starting from 0 at the root and increasing as nodes get topologically away from it. Interior nodes are arranged according to depth which is color coded from zero to 6 (horizontal legend). Terminals are color coded according to their clockwise order (vertical legend). The yellow boxes mark subtrees corresponding to equivalent parts in different objects.

Figure 14 :

 14 Figure 14: Subtrees that are geometrically accepted for wing node generation.

Figure 15 :

 15 Figure 15:The minimum bounding boxes of terminal nodes. Note that terminals A and B have similar bounding boxes, but the eccentricity of A is less than 1 whereas that of B is closer to 3.

Figure 16 :

 16 Figure 16: The pockets of void areas between a shape and its convex hull (shaded with color).

Figure 17 :

 17 Figure 17: Sample retrieval results on Kimia-216. The top row is the query shape. Each match is illustrated by the query object and the retrieved one where pairs of corresponding terminal segments are shaded by the same color.

  0.5 angle (in degrees) of iDM 0.1

Figure 18 :

 18 Figure 18: The retrieval results of objects selected from Kimia-99 with occluded parts. Each match is illustrated by the query object and the retrieved one where pairs of corresponding terminal segments are shaded by the same color.

Table 1 :

 1 The parameters' and distances' weights.

	ϵ	1%	th large 20% k long 10
	µ 0 0.75 th bent	7%	c bent 0.9
	µ 1 0.9

Table 2 :

 2 Threshold values

Table 3 :

 3 Terminals (N T ) Wings (N W ) Start points (N start ) Mean values of the factors that affect the complexity of our algorithm. The first two rows are for Kimia-99 and rows 3 and 4 are for Kimia-216 dataset.

	Wings of Terminal
	(N W oT )

Table 4 :

 4 The precision table for Kimia-99 dataset. It is made of 9 classes with 11 objects each. The total retrieval result is 98.18% in the nearest 10 matches. The highest two retrieval results are highlight in boldface.

	Method	1 st	2 nd 3 rd	4 th	5 th	6 th	7 th	8 th	9 th	10 th Total
	Shape context [6]	97	91	88	85	84	77	75	66	56	37	76.36
	Shock graph [33]	99	99	99	98	98	97	96	95	93	82	96.56
	Path similarity [4]	99	99	99	99	96	97	95	93	89	73	94.84
	Symbolic representation [11]	99	99	99	98	99	98	98	95	96	94	98.48
	MGVF [17]	99	99	99	99	99	99	99	98	90	87	97.77
	Two strategies [38]	99	99	99	98	99	99	99	97	96	84	97.87
	Height functions [42]	99	99	99	99	98	99	99	96	95	88	98.08
	Our method with the original CAT	98	95	95	90	89	84	78	67	66	50	82.02
	Our method with the extended CAT	99	99	99	99	99	99	99	98	96	85	98.18

Data: Shape A with n terminals and Shape B with m terminals Result: A pairing between the visual parts of the shapes and a cost estimate for the match. minCost = infinity; forall the terminal nodes T in A do if T NOT valid start point (Algorithm 3) then continue; end rotate the list of terminals so that T is the first; Algorithm 4: setup an n × m matrix of Options ; for i = n; i > 0; i-do for j = m; j > 0; j-do Algorithm 5: select the minCostOption at entry (i, j); end end if minCost > minCost at entry [START_REF] Cagri Aslan | Disconnected skeleton: Shape at its absolute scale[END_REF][START_REF] Cagri Aslan | Disconnected skeleton: Shape at its absolute scale[END_REF] then minCost = minCost at entry (1,1); end end /*Now minCost contains the estimated cost of the optimal correspondence between shapes A and B*/ return minCost; Algorithm 2: A rough pseudocode describing the dynamic time warping method that finds a minimal cost correspondence between the visual parts of two shapes A and B. Algorithm 3: Pseudocode that determines if a terminal is a valid start point. Terminals are always ordered in their anti-clock direction appearance along the boundary of the shape.

Data:

Data: Shape A with n terminals and Shape B with m terminals

Result: A matrix of options where every entry contains the options of matching terminal to terminal, terminal to wing, wing to terminal and wing to wing. forall the entry i, j in the matrix do add option pair: terminal i of A and terminal j of B; forall the wing nodes w of terminal i do add option pair: wing w of A and terminal j of B; end forall the wing nodes w of terminal j do add option pair: wing w of B and terminal i of A; end forall the wing nodes w1 of terminal i do forall the wing nodes w2 of terminal j do add option pair: wing w1 of A and wing w2 of B; end end //options to discard each terminal: add option: penalty(terminal i in A); add option: penalty(terminal j in B); end /*An "Option" contains information about matching two nodes (node A Option and node B Option )from the two shapes. It also saves information about the disjoint previous and next entries in the matrix. This information is essential to avoid having intersecting wings in the resultant correspondence*/ Algorithm 4: The pseudocode that sets up the matrix of options Data: A list of options at entry (i, j). Result: The minimal cost option "minCostOption" at this entry. This code also collects a path of pairs that constitutes the correspondence in the block [(i, j), (n, m)] of the matrix of options where (n×m) are the dimension sizes of the matrix. minCostOption = -1; minCost = infinity; forall the Option O at entry i, j do // the pair of nodes of O are node A O from shape A and node B O from shape B; fetch Prev-i and Prev-j;//subsequent disjoint entry in the matrix subOption = minCostOption in (Prev-i, Prev-j);

∪ path to subOption; end end /*The set of pairs of the minimal cost correspondence are found in the path to minCostOption at entry (1,1) in the matrix of options.*/ Algorithm 5: A rough pseudocode describing the optimal option selection at entry (i, j) of the matrix of options. Details related to the code that prevents intersections between wings that are selected into the correspondence are dropped to make the code more readable. rather simple building blocks of shape representation and matching, our approach produces good retrieval results on benchmark databases.

Having assessed the validity of our method on 2D image datasets, our future direction is to apply it on sketch-based 3D object retrieval. Generalizing this part correspondence matching to 3D [START_REF] Van Kaick | A survey on shape correspondence[END_REF] is an interesting problem. A possible direction for future work is to build a classifier based on the codebook of visual parts learnt from matching objects of the same class. A similar idea has been proposed by Kim et al. [START_REF] Vladimir | Learning part-based templates from large collections of 3d shapes[END_REF] in the context of 3D objects.

Acknowledgements

Zahraa Yasseen is supported by a scholarship from the Lebanese National Council for Scientific Research (CNRS). This work is also supported by a research grant from CNRS awarded to Ahmad Nasri.