Évelyne Touchette

Aude Chollet

Cédric Galéra

Eric Fombonne

Bruno Falissard

Michel Boivin

Maria Melchior

PhD Evelyne Touchette
email: eveline.touchette.1@ulaval.ca

Prior sleep problems predict internalising problems later in life

Keywords: Sleep problems, depressive and anxiety symptoms, longitudinal study, epidemiology

Background: One possible risk marker of later internalising symptoms is poor sleep, which is a problem for up to 40% of children. The present study investigated whether prior sleep problems could predict internalising symptoms over a period of 18 years of follow-up.

Methods:

The study sample included 1,503 French young adults from the TEMPO cohort (mean age=28.8±3.6 years) whose parents participate in the GAZEL cohort study. All TEMPO participants previously took part in a study of children's mental health and behaviour in 1991 (mean age=10.3±3.6 years) and 1999 (mean age=18.8±3.6 years). Sleep problems and internalising symptoms (depression, anxiety, somatic complaints) were assessed three times (1991, 1999, 2009) using the Achenbach System of Empirically Based Assessment (ASEBA) questionnaire. The association between sleep problems in 1991 and trajectories of internalising problems from 1991 to 2009 was tested in a multinomial logistic regression framework, controlling for sex, age, baseline temperament, behavioural problems and stressful life events, as well as family income, and parental history of depression.

Results:

We identified four trajectories of internalising symptoms: high-persistent (2.5%), high-decreasing (11.4%), low-increasing (11.6%), and low-persistent (74.5%). After controlling for covariates, compared to participants who did not have sleep problems in 1991, those who did were 4.51 times (95% CI = 1.54 -13.19, P=.006) more likely to have highpersistent internalising symptoms and 3.69 times (95% CI = 2.00 -6.82, P<.001) more likely to have high-decreasing internalising symptoms over the 18-year follow-up.

Limitations: Sleep problems and internalising symptoms were based on self-report questions, results should be interpreted with due caution.

Conclusions: Sleep problems early in life are associated with an increased likelihood of internalising symptoms that persist from childhood to adulthood.

Introduction

Internalising problems, defined as a spectrum of symptoms of depression, anxiety, withdrawal, and somatic complaints [START_REF] Achenbach | Replication of empirically derived syndromes as a basis for taxonomy of child/adolescent psychopathology[END_REF], affect up to 40% of youth [START_REF] Crawford | Internalizing symptoms in adolescents: Gender differences in vulnerability to parental distress and discord[END_REF]. In up to 50% of cases, these symptoms persist over time, affecting quality of life, interpersonal relationships, as well as academic and occupational attainment [START_REF] Bongers | The normative development of child and adolescent problem behavior[END_REF]. It is therefore important to identify youth at risk of persistent internalising symptoms in order to address their mental health needs early on and potentially prevent a long-term cycle of poor outcomes.

One possible risk marker of later internalising symptoms is poor sleep, which is a problem for up to 40% of children [START_REF] Mindell | Developmental features of sleep[END_REF]. Although sleep problems are sometimes considered a sign of anxiety and depression [START_REF] Ostberg | Feeding and sleeping problems in infancy--a follow-up at early school age[END_REF], research shows that they often precede the occurrence of other manifestations of psychological difficulties (such as excessive worry, unexplained somatic symptoms, sadness or social withdrawal). Sleep problems are associated with increased risk of internalising symptoms (particularly anxiety) in toddlers [START_REF] Jansen | Does disturbed sleeping precede symptoms of anxiety or depression in toddlers? The generation R study[END_REF], children [START_REF] Pesonen | Sleep duration and regularity are associated with behavioral problems in 8-year-old children[END_REF], adolescents [START_REF] Gregory | Sleep problems in childhood: a longitudinal study of developmental change and association with behavioral problems[END_REF], and adults [START_REF] Gregory | Prospective longitudinal associations between persistent sleep problems in childhood and anxiety and depression disorders in adulthood[END_REF]. For instance, a longitudinal study which followed 4-19 year olds over a period of 14 years, found that various types of sleep problems, as reported by the parents using the Child Behavioural Checklist (CBCL) (e.g., "decreased sleep", "overtiredness", and "trouble sleeping") predicted later depression and anxiety [START_REF] Gregory | Parent-reported sleep problems during development and self-reported anxiety/depression, attention problems, and aggressive behavior later in life[END_REF]. However, to our knowledge, the relationship between sleep problems early on and trajectories of psychological difficulties has not been studied. In the present study, we examined the relationship between sleep problems measured between 4 and 16 years of age, and internalising symptoms over an 18 year follow-up period, accounting for factors which could explain this association (sex [START_REF] Hankin | Sex differences in adolescent depression: stress exposure and reactivity models[END_REF], temperament [START_REF] Touchette | Factors associated with fragmented sleep at night across early childhood[END_REF], behavioural difficulties [START_REF] Silk | Adolescents' emotion regulation in daily life: links to depressive symptoms and problem behavior[END_REF], stressful life events [START_REF] Eley | Specific life events and chronic experiences differentially associated with depression and anxiety in young twins[END_REF], family socioeconomic position [START_REF] Arber | Gender and socio-economic patterning of self-reported sleep problems in Britain[END_REF], [START_REF] Lorant | Socioeconomic differences in psychiatric in-patient care[END_REF][START_REF] Melchior | Socioeconomic position predicts long-term depression trajectory: a 13-year follow-up of the GAZEL cohort study[END_REF], and parental depression [START_REF] Swanson | Maternal mood and sleep in children of women at risk for perinatal depression[END_REF].

Methods

Study population

The TEMPO (Trajectoires Epidémiologiques en Population) cohort study, based in France, was set up in 2009 to follow-up young adults (22-35 years) who had taken part in a study of children's psychological problems and access to mental health care in 1991. The original sample of children surveyed in 1991 (GAZEL Youth study, n=2,498) was selected among 4-16 year olds whose parents participate in the GAZEL cohort study which follows 20,624 employees of a large French public-sector utility company by yearly mailed questionnaire since 1989 [START_REF] Goldberg | Cohort profile: the GAZEL Cohort Study[END_REF]. The 1991 sample of children was selected to match key characteristics of children in France (number of children per family and occupational grade of head of household) (Fombonne & Vermeersch, 1997a, 1997b). Data on participating children were collected via parental reports. In 1999, data were collected by parental reports (n=1,268) and youth self-reports (n=1,148).

In 2009, all living parents of children who took part in the GAZEL Youth Study in 1991 received a letter asking them to forward the TEMPO study questionnaire to their son/daughter. Between 1991 and 2009, 16 participants died and 4 were too ill or disabled to answer. The overall response rate to the TEMPO questionnaire was 44.5% (n=1,103), which is comparable to response rates of other mental health surveys in France [START_REF] Alonso | Sampling and methods of the European Study of the Epidemiology of Mental Disorders (ESEMeD) project[END_REF].

Leading reasons for non-participation were non-transmission of the questionnaire by the parent (34.4%) or the youth's lack of interest (28.5%). Compared to 2009 respondents, non-respondents were older, more likely to have parents who were divorced, and had low socioeconomic backgroundbut did not vary with regard to their parents or their own overall psychological characteristics. The TEMPO study was reviewed and approved by an ethical board (CCTIRS: Comité Consultatif sur le Traitement des Informations pour la Recherche en Santé) and France's national committee for data protection (CNIL: Commission Nationale Informatique et Liberté).

Outcome measures

Participants' internalising symptoms were measured three times (1991,1999,2009) using the Achenbach System of Empirically Based Assessment (ASEBA). This widely used instrument comprises 118 items which assess behaviour (internalising and externalising symptoms) over a six-month period [START_REF] Achenbach | Manual for the Teacher's report form and 1991 profile[END_REF]. This questionnaire has previously been validated in France [START_REF] Fombonne | The use of questionnaires in child psychiatry research: measuring their performance and choosing an optimal cut-off[END_REF]; [START_REF] Stanger | Epidemiological comparisons of American and French children: Parent reports of problems and competencies for ages 6-11[END_REF]. As suggested by Gregory et al. [START_REF] Gregory | Sleep items in the child behavior checklist: a comparison with sleep diaries, actigraphy, and polysomnography[END_REF], we used the most rigorous CBCL sleep item: "Does your child have sleep problems?" which has been shown to correlate with sleep latency assessed both by diary (P=.008) and actigraphic measures (P=.029).We considered sleep problems to be present when the parents reported them "sometimes" or "often" (vs. never). In additional analyses, we used a stricter cut-off, comparing individuals who had sleep problems "often" to those who had sleep problems "sometimes" or "never". Internalising symptoms (i.e. depressive symptoms, anxiety, withdrawal, and somatic complaints) were assessed by parental reports (1991: n=32 items), by parental and youth self-reports combined (1999: n=31 items) and by young adult self-reports (2009: n=43 items).Each time, we excluded "sleep problems" from the items included to construct the internalising symptoms score.

Internalising symptoms scores had high internal consistency at all three measurement points: Cronbach's alphas as follows: 1991=0.83, 1999=0.88, and 2009=0.93. To estimate trajectories of internalising symptoms, we used continuous scores. Additionally, we also dichotomised internalising symptoms to identify participants with clinically significant symptoms levels, using the 85 th percentile score suggested by Amone-P 'Olak et al.(Amone-P'Olak et al., 2009).

Description of covariates

In the present study, we controlled for the following covariates:

1) Sex (male vs. female);

2) Age (studied as a continuous variable);

3) Childhood unstable temperament, reported retrospectively for age 7-10 years by parents at time 2using the French version of the Emotionality Activity Sociability (EAS) questionnaire [START_REF] Gasman | Cross-cultural assessment of childhood temperament. A confirmatory factor analysis of the French Emotionality Activity and Sociability (EAS) questionnaire[END_REF]. Unstable temperament score was defined based on answers to five items: "Cried easily?", "Was too sensitive?", "Got excited over nothing?", "Was easily upset?", and "Reacted strongly when upset?". This score was then standardized to a mean of 0 and standard deviation of 1. 4) Initial externalizing problems (studied as a continuous variable) comprising aggressive and rule-breaking behaviours of the ASEBA questionnaire in 1991; 5) Stressful life events were measured by a combination of experiences and situations reported in 1991 (i.e., school difficulties, parental stress, childhood illness, childhood social isolation, illness of a close family member/friend, move, parental quarrels, death of a close family member/friend, parental unemployment/financial problems, parental absence from home) and parental divorce prior to age 18 as assessed based on parents' reports in the yearly GAZEL study questionnaire. We summed all stressful life events and parental divorce and studied a total score ranging from 0 to 11; 6) Family socio-economic status was assessed by parent-reported family income in 1989 (dichotomized at the median value of 1,981 €/month) (INSEE, 2009) 7) Parental depression was assessed using two sources of information: 1) ≥2 self-reports of depression in the yearly GAZEL study questionnaire , and 2) TEMPO participants' 2009 reports of their parents' lifetime experience of depression (yes vs. no) on the National Institute of Health-Family Inventory for Genetic Studies (NIH-FIGS)

questionnaire [START_REF] Maxwell | Manual for the Family Interview for Genetics Studies (FIGS)[END_REF].

Statistical Analyses

The analyses were based on 1,503 TEMPO participants who took part in at least two out of three study assessments (1991,1999,2009). First, we tested age and sex-adjusted concomitant associations between sleep problems and clinically significant internalising difficulties. Logistic regression models were used to test whether sleep problems in 1991 predicted the emergence of later internalising difficulties after adjusting for internalising difficulties in 1991 and the aforementioned covariates.

Next, we identified trajectories of internalising difficulties from 1991 to 2009 using a semiparametric modeling method [START_REF] Nagin | Group-based modeling of development[END_REF] implemented in SAS® (SAS Institute Inc., Cary, NC, PROC TRAJ). Based on the maximum Bayesian information criterion (BIC) the model identified four trajectories as the best fit to the data. Each individual was assigned to a specific trajectory based on a posterior probability of belonging to that group. We tested univariate associations between both sleep problems in 1991 and trajectories of internalising symptoms and all potential covariates using chi-square tests (for categorical variables), and ttests or ANOVAs (for continuous variables).

Finally, a multinomial multivariate regression was performed to assess the association between prior sleep problems and trajectories of internalising symptoms adjusting for participants' age, sex, and all covariates associated with the presence of sleep problems in 1991 and/or trajectories of internalising symptoms (P<0.10). Analyses were performed using the SPSS statistical software (version 16.0, SPSS Inc, Chicago, ILL).

Results

Transversal analyses between sleep and internalising problems

In the present study, 17.4% (n=260) of children had sleep problems in 1991, 29.9% (n=330) in 1999 and 34.3% (n=364) in 2009. Restricting the definition to frequent sleep problems, prevalences were as follows: 3.7% (n=56) in 1991, 7.1% (n=78) in 1999 and 7.5% (n=80) in 2009. Table 1 shows moderate age and sex-adjusted correlations between sleep problems in 1991 and internalising problems in 1991 (r=.18), 1999 (r=.17), and no association with internalising symptoms in 2009 (r=.03). As shown in Table 2, adjusting for internalising problems in 1991 and all covariates, sleep problems in 1991 were associated with internalising problems in 1999 (OR=2.18, 95% CI=1.22 -3.89, P=.009) but not in 2009 (OR=0.71, 95% CI=0.38-1.39, P=.33) [Insert Tables 1 and2 about here] had high internalising symptoms in 1991 which decreased during follow-up (high-decreasing trajectory). In addition, 11.6% (n=176) had low levels of internalising symptoms in 1991 which increased over the course of time (low-increasing trajectory). Finally, 2.5% of study participants (n=37) had high levels of internalising symptoms throughout follow-up (highpersistent trajectory).

Trajectories of internalising symptoms

[Insert Figure 1 about here]

Selection of covariates

As shown in Table 3, sleep problems in 1991 were significantly associated with stressful life events (P<.001), externalizing problems (P<.001), an unstable temperament (P=.001) and parental depression (P=.01). Trajectories of internalising symptoms were significantly associated with participants' age (P<.001), sex (P<.001), stressful life events (P<.001), externalizing problems (P<.001), an unstable temperament (P<.001), and parental depression (P<.001), but not with family income (P=.20). For consistency with prior studies, family income was nevertheless controlled for.

[Insert Table 3 about here]

Sleep problems in 1991 and trajectories of internalising symptoms

We found an association between sleep problems in 1991 and trajectories of internalising symptoms (P<.001) and this association was statistically significant in female and male participants (respectively P<0.001 and P<0.001). Sleep problems in 1991 were more frequent in participants in the high-persistent (girls=44.8%, boys=57.1%) or high-decreasing (girls=41.4%, boys=38.6%) trajectories of internalising symptoms than in the low-increasing (girls=17.5%, boys=18.0%) and low-persistent (girls=14.0%, 11.8%) trajectories.

Associations between sleep problems and internalising symptoms trajectories did not vary with sex (interaction term not statistically significant).

Multivariate multinomial regression models testing associations between sleep

problems in 1991 and trajectories of internalising symptoms adjusting for all covariates are presented in Table 4. As expected, compared to participants with no sleep problems, those who had sleep problems in 1991 were more likely to have a high-persistent trajectory of internalising symptoms (OR=4.51, 95% CI: 1.54-13.19, P=0.006) or a high-decreasing trajectory (OR=3.69, 95% CI = 2.00 -6.82, P<.001). We found no association between sleep problems in 1991 and a low-increasing internalising symptoms trajectory (P=.89). Additional variables associated with the high-persistent trajectory of internalising symptoms were female sex and an unstable temperament. Additional variables associated with the high-decreasing internalising symptoms trajectory were female sex, an unstable temperament, externalizing problems, and stressful life events [Insert Figure 2 and Table 4 about here]

Discussion

Studying a sample of young adults drawn from the community, we found that participants who had prior sleep problems had an elevated likelihood of experiencing longterm internalising problems. After accounting for important covariates, children who had sleep problems were 4.5 times more likely to have persistently high levels of internalising symptoms into young adulthood but, 3.7 times more likely to have high levels of internalising symptoms that later decreased. This association did not appear to be sex-specific; it was equivalently observed in girls and boys.

The results of the present study corroborate the predictive link between prior sleep problems and internalising problems later in life observed in previous longitudinal studies [START_REF] Gregory | Sleep problems in childhood: a longitudinal study of developmental change and association with behavioral problems[END_REF][START_REF] Jansen | Does disturbed sleeping precede symptoms of anxiety or depression in toddlers? The generation R study[END_REF]. Dyssomnias, parasomnias, and short sleep duration at 18 months of age increase the risk of occurrence of depressive or anxiety symptoms at age 3 even after taking into account pre-existing depressive or anxiety symptoms [START_REF] Jansen | Does disturbed sleeping precede symptoms of anxiety or depression in toddlers? The generation R study[END_REF]. Moreover, sleep problems at age 4 predict behavioral/emotional problems in mid-adolescence [START_REF] Gregory | Sleep problems in childhood: a longitudinal study of developmental change and association with behavioral problems[END_REF]. However, we found no association between early sleep problems in 1991 and internalising problems in young adulthood, implying that internalising symptoms in adulthood have other determinants. In order to clarify the mechanisms underlying this association, future studies should focus on prior sleep problems measured objectively and mental health assessed later in life.

ASEBA measures symptoms present over the last 6 months [START_REF] Achenbach | Manual for the Teacher's report form and 1991 profile[END_REF]. This suggests that the ASEBA may reflect trait as well as state expressions of psychopathology [START_REF] Biederman | CBCL clinical scales discriminate ADHD youth with structured-interview derived diagnosis of oppositional defiant disorder (ODD)[END_REF]. For example, individuals in the high-persistent internalising symptoms trajectory could indicate more a trait than a state of internalising symptoms compared with the other internalising symptoms trajectories. Different mechanisms could explain why prior sleep problems predict youths' internalising problems. First, this association could reflect common genetic influences, as suggested by a study of three hundred 8-year-old twin pairs [START_REF] Gregory | Associations between sleep problems, anxiety, and depression in twins at 8 years of age[END_REF]. Several physiological pathways through the monoaminergic system (serotonin, norepinephrine, and dopamine) and potentially the glutamatergic system [START_REF] Hashimoto | Emerging role of glutamate in the pathophysiology of major depressive disorder[END_REF], [START_REF] Willner | Dopamine and depression: a review of recent evidence. I. Empirical studies[END_REF] are involved in both sleep and emotional problems.. In addition, chronic stress induced by persistent sleep problems could affect the hypothalamic-pituitary-adrenal (HPA) axis [START_REF] Buckley | On the interactions of the hypothalamic-pituitaryadrenal (HPA) axis and sleep: normal HPA axis activity and circadian rhythm, exemplary sleep disorders[END_REF] and result in an unusual pattern of cortisol secretion [START_REF] Feder | Twenty-four-hour cortisol secretion patterns in prepubertal children with anxiety or depressive disorders[END_REF], which in turn, may have an effect on higher-order regulatory systems involved in emotional regulation. Similar HPA axis deregulations have been found in humans with depression and in prenatally stressed animals, suggesting that gestational stress could increase the risk of developing internalising problems [START_REF] Weinstock | Does prenatal stress impair coping and regulation of hypothalamicpituitary-adrenal axis?[END_REF]. Moreover, sleep problems decrease daytime vigilance [START_REF] Peters | The sensitivity of a PDA-based psychomotor vigilance task to sleep restriction in 10-yearold girls[END_REF] which is increasingly thought to be associated with disproportionate levels of behavioral/emotional problems [START_REF] Lavigne | Sleep and behavior problems among preschoolers[END_REF][START_REF] Stein | Sleep and behavior problems in school-aged children[END_REF].

Conversely, shared environmental factors seem to explain a moderate proportion of the variance in the association between sleep problems and psychological problems such as depression [START_REF] Gregory | Associations between sleep problems, anxiety, and depression in twins at 8 years of age[END_REF]. Particularly, inadequate parental behaviors impeding selfsoothing abilities could be implicated into both sleep-wake and emotional regulatory systems [START_REF] Touchette | Risk factors and consequences of early childhood dyssomnias: New perspectives[END_REF]. Underlying child characteristics, such as temperament, might also be involved. In the present study, unstable childhood temperament was independently associated with later internalising symptoms but did not explain the association between prior sleep problems and the emergence of internalising symptoms 8 years later.

Several methodological strengths and limitations should be mentioned. We studied a longitudinal sample of community-based children for a period of 18 years where prior sleep problems were measured independently from internalising symptoms and analyses were adjusted for pre-existing individual and familial characteristics. As in many epidemiological studies, our measure of sleep was based on a single CBCL sleep item, which is well correlated with sleep latency assessed by diary and actigraphic measures [START_REF] Gregory | Sleep items in the child behavior checklist: a comparison with sleep diaries, actigraphy, and polysomnography[END_REF]. As others [START_REF] Gregory | Prospective longitudinal associations between persistent sleep problems in childhood and anxiety and depression disorders in adulthood[END_REF][START_REF] Gregory | Parent-reported sleep problems during development and self-reported anxiety/depression, attention problems, and aggressive behavior later in life[END_REF], we used parental perceptions of sleep problems, which could be different from polysomnographic sleep problems. In addition, we used ASEBA scores to assess psychological difficulties, which are less specific than psychiatric diagnoses. Nevertheless, the ASEBA system has good psychometric properties [START_REF] Carter | Assessment of young children's socialemotional development and psychopathology: recent advances and recommendations for practice[END_REF] and has frequently been used to screen for clinically significant depressive/anxiety problems throughout the life course.

Conclusion

Sleep problems early in life may contribute to the development of internalising problems which persist for a significant period of life. Specifically, children with sleep problems appear vulnerable to internalising symptoms in adolescence, which sometimes persist into young adulthood. These results highlight the importance of investigating the presence of sleep problems in children, as these may be indicative of risk of long-term psychological difficulties. Future studies should explore the biological mechanisms underlying the association between sleep problems and internalising symptoms and test whether treatment of sleep problems can help contain the risk of internalising psychological difficulties later on in life.

Figure 1

 1 Figure 1 depicts trajectories of internalising symptoms between 1991 and 2009. The


Julien Tierny
email: tierny@lifl.fr

Jean-Philippe Vandeborre
email: vandeborre@lifl.fr

Mohamed Daoudi
email: daoudi@lifl.fr

3D Mesh Skeleton Extraction Using Topological and Geometrical Analyses

Mapping function. (c) Dual Reeb graph. (d) Constriction approximation. (e) Constriction enhanced graph. (f) Application to deformation.

Without pre-processing stages and without input parameters, our method provides nice-looking and affineinvariant skeletons, with satisfactory execution times. This makes enhanced topological skeletons good candidates for applications needing high level shape representations, such as mesh deformation (experimented in this paper), retrieval, compression, metamorphosis, etc.

Abstract

This paper describes a novel and unified approach for Reeb graph construction and simplification as well as constriction approximation on 3D polygonal meshes. The key idea of our algorithm is that discrete contours -curvescarried by the edges of the mesh and approximating the continuous contours of a mapping function -encode both topological and geometrical shape characteristics.

Firstly, mesh feature points are computed. Then they are used as geodesic origins for the computation of an invariant mapping function that reveals the shape most significant features. Secondly, for each vertex in the mesh, its discrete contour is computed. As the set of discrete contours recovers the whole surface, each of them can be analyzed, both to detect topological changes or constrictions. Constriction approximations enable Reeb graphs refinement into more visually meaningful skeletons, that we refer as enhanced topological skeletons.

Introduction

Polygonal mesh is a widely used representation of 3D shapes, mainly for exchange and display purposes. However, many applications in computer graphics need higher level shape descriptions as an input. Topological skeletons have shown to be interesting shape descriptions [START_REF] Biasotti | An overview on properties and efficacy of topological skeletons in shape modelling[END_REF]. They benefit diverse fields like shape metamorphosis, deformation [START_REF] Lewis | Pose space deformations: A unified approach to shape[END_REF], retrieval [START_REF] Hilaga | Topology matching for fully automatic similarity estimation of 3D shapes[END_REF], texture mapping, etc.

Many topological approaches study the properties of real valued functions computed over triangulated surfaces. Most of the time, those functions are provided by the application context, such as scientific data analysis [START_REF] Bremer | Topological hierarchy for functions on triangulated surfaces[END_REF]. When dealing with topological skeletons, it is necessary to define an invariant and visually interesting mapping function, which remains an open issue [START_REF] Biasotti | An overview on properties and efficacy of topological skeletons in shape modelling[END_REF].

Moreover, traditional topological graph construction algorithms assume that all the information brought by the mapping function is pertinent, while in practice, this can lead to large graphs [START_REF] Ni | Fair Morse functions for extracting the topological structure of a surface mesh[END_REF][START_REF] Carr | Simplifying flexible isosurfaces using local geometric measures[END_REF], encoding noisy details.

Finally, topological approaches cannot discriminate visually interesting sub-parts of identified connected components, like the phalanxes of a finger. This is detrimental to certain applications, such as mesh deformation.

In this paper, we propose a novel and unified method which addresses the above issues. Given a connected triangulated surface T ,featurepointsarefirstlyextracted(fig. 1(a)) in order to compute an invariant mapping function, noted f m (fig. 1(b)), which reveals the shape most signifi-cant parts. Secondly, for each vertex in the mesh, we compute its discrete contour,aconnectedcurvetraversingitand locally minimizing f m gradient. We show that a topological analysis of those discrete contours enables a pertinent Reeb graph construction and simplification (fig. 1(c)), without any input parameter. Finally, we show that a geometrical analysis of discrete contours can approximate constrictions on prominent components (fig. 1(d)), enabling the refinement of Reeb graphs into enhanced topological skeletons (fig. 1(e)).

This paper is structured as follows. Firstly, we introduce topological skeleton related work. Secondly, we define our mapping function f m .Thirdly ,wepresentouralgorithmfor discrete contour computation, which is used both for the Reeb graph construction and simplification as well as the constriction approximation. Finally, we comment on experimental results and discuss about possible applications, like mesh deformation (fig. 1(f)).

Related work

Several approaches have been explored for the decomposition of polygonal meshes into meaningful subcomponents, to extract skeletal representations of shapes.

In comparison to mesh segmentation [START_REF] Shamir | A formalization of boundary mesh segmentation[END_REF] and traditional skeleton extraction [START_REF] Blum | Shape description using weighted symmetric axis features[END_REF][START_REF] Wu | Domain connected graph: the skeleton of a closed 3D shape for animation[END_REF], topological approaches, based on Morse and Reeb graph theories [START_REF] Morse | Relations between the critical points of a real function of n independant variables[END_REF][START_REF] Reeb | Sur les points singuliers d'une forme de Pfaff complètement intégrable ou d'une fonction numérique[END_REF][START_REF] Milnor | [END_REF], present the advantage to preserve the topological properties of the shape [START_REF] Biasotti | An overview on properties and efficacy of topological skeletons in shape modelling[END_REF] (number of loops, number and relations between components, etc.). However, with regard to shape skeletons, we identify three main drawbacks in topological approaches, successively addressed in this paper.

Firstly, it is difficult to define an invariant and visually interesting mapping function. Secondly, constructing and transforming a topological graph into a manageable skeleton is not a trivial problem. Finally, topological approaches decompose a surface into connected sub-components only. This means that visually interesting sub-parts of identified connected components will not be discriminated: for example, a finger of a hand model will not be decomposed into phalanxes.

Mapping functions

Differential topology based approaches study the properties of real valued functions, that we refer as mapping functions,d e fi n e do ni n p u ts u r f a c e s ,e i t h e rt oc o n s t r u c tR e e b graphs [START_REF] Shinagawa | Surface coding based on morse theory[END_REF][START_REF] Cole-Mclaughlin | Loops in Reeb graphs of 2-manifolds[END_REF], contour trees [START_REF] Carr | Simplifying flexible isosurfaces using local geometric measures[END_REF], level set diagrams [START_REF] Lazarus | Level set diagrams of polyhedral objects[END_REF] or Morse complexes [START_REF] Ni | Fair Morse functions for extracting the topological structure of a surface mesh[END_REF][START_REF] Bremer | Topological hierarchy for functions on triangulated surfaces[END_REF]. Those functions are often brought by the application context: terrain modeling [START_REF] Shinagawa | Surface coding based on morse theory[END_REF], MRI analysis [START_REF] Carr | Simplifying flexible isosurfaces using local geometric measures[END_REF], molecular analysis [START_REF] Bremer | Topological hierarchy for functions on triangulated surfaces[END_REF], etc.

When dealing with topological skeletons, it is necessary to define a scalar function which satisfies invariance and stability constraints, and which also affords a topological description that highlights visually significant surface subcomponents.

Lazarus and Verroust [START_REF] Lazarus | Level set diagrams of polyhedral objects[END_REF] introduced such a function, defined by the geodesic distance (the length of the shortest path between vertices) from a source vertex to any other vertex in the mesh. It leads to visually interesting results for natural objects because it is invariant to geometrical transformations and it is robust against variations in model pose [START_REF] Katz | Mesh segmentation using feature point and core extraction[END_REF]. Due to a lack of stability, within the framework of shape retrieval, Hilaga et al. [START_REF] Hilaga | Topology matching for fully automatic similarity estimation of 3D shapes[END_REF] proposed to integrate this function all over the mesh. Unfortunately, from our experience, that function generates an important amount of critical points, configurations where the gradient of the function vanishes, which makes the construction of visually meaningful graphs more complex.

In our method, to reveal the shape most significant features, we focus on feature points. Feature points are mesh vertices located on extremities of prominent components [START_REF] Katz | Mesh segmentation using feature point and core extraction[END_REF]. Mortara and Patanè [START_REF] Mortara | Affine-invariant skeleton of 3D shapes[END_REF] proposed to select as feature points the vertices where Gaussian curvature exceeds a given threshold, but this cannot resolve extraction on constant curvature areas. Katz et al. [START_REF] Katz | Mesh segmentation using feature point and core extraction[END_REF] developed an algorithm based on multi-dimensional scaling in quadratic execution complexity. In this paper, we propose a robust and straightforward algorithm for feature point extraction (fig. 1(a)). Moreover, we use them as geodesic origins for the definition of our mapping function (fig. 1(b)). Such a function well reveals the most visually significant parts of the mesh, generating manageable critical point sets.

Graph construction and simplification

AReebgraph [START_REF] Reeb | Sur les points singuliers d'une forme de Pfaff complètement intégrable ou d'une fonction numérique[END_REF]isatopologicalstructurethatencodes the connectivity relations of the critical points of a scalar function defined on an input surface. More formally, Reeb graphs are defined as follows:

Definition 1 (Reeb graph) Let f : M → R be a scalar function defined on a compact manifold M .TheReebgr aph of f is the quotient space of f in M × R by the equivalence relation (p 1 ,f(p 1 )) ∼ (p 2 ,f(p2)),ifandonlyif:    f (p 1 )=f (p 2 ) p 1 and p 2 belong to the same connected component of f -1 (f (p 1 ))
Figure 2 gives an example of a Reeb graph computed on a bi-torus with regard to the height function and well illustrates the fact that Reeb graphs can be used as skeletons.

Constructing a Reeb graph from a scalar function f computed on a triangulated surface first requires to identify the set of vertices corresponding to critical points. With this aim, several formulations have been proposed [START_REF] Edelsbrunner | Simulation of simplicity: at e c h n i q u et oc o p ew i t hd e g e n e r a t ec a s e si ng e o m e t r i ca lgorithms[END_REF][START_REF] Takahashi | Algorithms for extracting correct critical points and constructing topological graphs from discrete geographical elevation data[END_REF] to identify local maxima, minima and saddles, observing for each vertex the evolution of f at its direct neighbors. Several algorithms have been developed to construct Reeb graphs from the connectivity relations of these critical points [START_REF] Cole-Mclaughlin | Loops in Reeb graphs of 2-manifolds[END_REF][START_REF] Carr | Simplifying flexible isosurfaces using local geometric measures[END_REF], most of them in O(n × log(n)) steps, with n the number of vertices in the mesh. However, they assume that all the information brought by the scalar function f is relevant [START_REF] Ni | Fair Morse functions for extracting the topological structure of a surface mesh[END_REF][START_REF] Carr | Simplifying flexible isosurfaces using local geometric measures[END_REF]. Consequently, they assume that all the identified critical points are meaningful, while in practice, this hypothesis can lead to unmanageably large Reeb graphs. To overcome this issue, Ni et al. [START_REF] Ni | Fair Morse functions for extracting the topological structure of a surface mesh[END_REF] developed a user-controlled simplification algorithm. Bremer et al. [START_REF] Bremer | Topological hierarchy for functions on triangulated surfaces[END_REF] proposed an interesting critical point cancellation technique based on a persistence threshold. Attene et al. [START_REF] Attene | Shape understanding by contour-driven retiling[END_REF] proposed as e d u c i n ga p p r o a c h ,u n i f y i n gt h eg r a p hc o n s t r u c t i o na n d simplification, but it is conditioned by a slicing parameter.

In this paper, we propose a discrete formulation of contours, connected subsets of level lines, which enables, without any input parameter, the construction of visually meaningful Reeb graphs (fig. 1(c)).

Constriction computation

Hétroy and Attali [START_REF] Hétroy | From a closed piecewise geodesic to a constriction on a closed triangulated surface[END_REF] define constrictions as simple closed curves, whose length is locally minimal. Constrictions are located on the narrowest parts of a surface. This notion is of major interest for segmenting individual subcomponents identified with standard topological approaches into more significant parts, for deformation purpose for example. Recently, Hétroy [START_REF] Hétroy | Constriction computation using surface curvature[END_REF] showed that constriction detection could be achieved by analyzing surface curvature.

In this paper, we propose to analyze the geometrical characteristics of discrete contours, and particularly their concavity, to approximate constrictions (fig. 1(d)), in order to decompose previously identified components into more visually interesting parts (fig. 1(e)).

Method overview

Given a connected triangulated surface T ,w ep r o p o s e in this paper a unified method to decompose T into visually meaningful sub-parts, considering the topological and geometrical characteristics of discrete contours.

The algorithm proceeds in three stages. Firstly, mesh feature points are extracted (fig. 1(a)) in order to compute an invariant and visually interesting mapping function (fig. 1(b)), denoted f m in the rest of the paper. Secondly, for each vertex in the mesh, we compute its discrete contour,a curve traversing it and approximating f m continuous contour. Finally, as the set of discrete contours recovers the entire mesh, it is possible to analyze each contour characteristics, either to detect topological changes (fig. 1(c)) or to detect curvature transitions (fig. 1(d)).

Our scientific contribution resides in three points. (i) We propose an original and straightforward algorithm for feature point extraction. It can resolve extraction on constant curvature areas -such as spheres -and it is robust against variations in mesh sampling and model pose. (ii) We show that a discrete contour formulation enables, without re-meshing and without any input parameter, a pertinent Reeb graph construction, providing visually meaningful graphs, affine-invariant and robust to variations in mesh sampling. (iii) We show that the geometrical information brought by discrete contours enables the approximation of constrictions on prominent components and consequently Reeb graph refinement.

Mapping function

To compute visually meaningful topological skeletons, we first have to define a mapping function that will highlight the most significant parts of the mesh. In order to fit application constraints, this function has to present stability and invariance properties. Geodesic distances are affineinvariant and robust to variations in model pose. From an algorithmic point of view, they can be approximated by the Moore-Dijkstra algorithm (distance minimizing in weighted graphs). In the rest of this paper, we will refer to δ(v 1 ,v 2 ) as the normalized approximation of the geodesic distance from vertex v 1 to vertex v 2 ,normalizedwithregard to mesh global extrema.

Feature point extraction

Feature points are mesh vertices located on extremities of prominent components. As they highlight the most significant features of the shape, they are used in our mapping function computation as origins for geodesic distance evaluation. To extract them, we propose a quite straightforward algorithm, based on topological tools. Most of the geodesic based mapping function local extrema appear at extremities of prominent components (figs. 3(a) and 3(b)), mainly because gradient vanishes in those configurations. Therefore, we propose to realize a crossed analysis, using two geodesic based mapping functions -whose origins are the mesh most distant vertices -and to intersect the sets of their local extrema. Let v s1 and v s2 be the most geodesic distant vertices of aconnectedtriangulatedsurfaceT ,computedwiththeTree Diameter algorithm [START_REF] Lazarus | Level set diagrams of polyhedral objects[END_REF]. In figure 3, v s1 is located at the extremity of the wrist (fig. 3(a)) while v s2 is located at the extremity of the middle finger (fig. 3(b)).

Let f g1 and f g2 be two scalar functions defined on each vertex v of T ,asfollows:

f g1 (v)=δ(v, v s1 ) (1) f g2 (v)=δ(v, v s2 )
(2) Basing on the critical point classification proposed in [START_REF] Cole-Mclaughlin | Loops in Reeb graphs of 2-manifolds[END_REF], a local minimum is defined as a vertex such that all its direct neighbors have an upper function value. Reciprocally, we define a local maximum as a vertex such that all its direct neighbors have a lower function value. Let E 1 be the set of local extrema (minima and maxima) of f g1 (in yellow in fig. 3(a)) and E 2 be the set of local extrema of f g2 (in cyan in fig. 3

F = E 1 ∩ E 2
(3) In practice, f g1 and f g2 local extrema which correspond to feature points do not appear exactly on the same vertices but in the same geodesic neighborhood.T h e r e f o r e ,t h ei ntersection constraint is relaxed as follows, with ∈ [0, 1] the radius of the geodesic neighborhood (geodesic distance are normalized):

v ∈ F ⇐⇒        ∃v e1 ∈ E 1 /δ (v, v e1 ) < ∃v e2 ∈ E 2 /δ (v, v e2 ) < δ(v, v fi ) > ∀v fi ∈ F ∈ [0, 1] (4)
From our experience, using only two geodesic mapping functions (f g1 and f g2 )a nds e t t i ng =0 .05 give accurate results. Moore-Dijkstra's algorithm is a time complexity bottleneck. f g1 and f g2 are computed each in O(n×log(n)) steps, with n the number of vertices in T .

In this paragraph, we presented a straightforward algorithm for mesh feature point extraction, in O(n × log(n)) steps, with n the number of vertices in T .Thisalgorithmis based on geodesic distance evaluations. Therefore it is stable and invariant to geometrical transformations and robust to variations in model pose. Furthermore, in order to select feature points, the mapping function gradient is analyzed.

No hypothesis is required about mesh sampling. Consequently, this algorithm is robust against variations in mesh sampling, as illustrated in figure 4: feature points are similar when the resolution decreases. Moreover, it achieves correct extraction on constant curvature areas, such as spheres, as shown in fig. 5(a).

Mapping function computation

When dealing with topological skeletons, it is necessary to define an invariant and visually interesting mapping function, which reveals the most significant parts of the mesh. Moreover, the mapping function should not generate an unmanageable set of critical points, in order to make the graph simplification easier. From our experience, this is not the case of the function presented in [START_REF] Hilaga | Topology matching for fully automatic similarity estimation of 3D shapes[END_REF].

Firstly, to guarantee invariance to geometrical transformations and robustness against variations in model pose, geodesic distances are used. Secondly, to define a visually interesting mapping function, feature points are taken as origins for geodesic distance evaluations. Therefore, we propose the following mapping function, noted f m in the rest of the paper, which computes in each vertex v of T the geodesic distance to the closest feature point:

f m (v)= f c (v) -min v∈T f c (v) max v∈T f c (v) -min v∈T f c (v) (5)
f m is a normalized version of the function f c ,definedas follows (f c (v) ≥ 0, ∀v ∈ T ):

f c (v)=1-δ ′ (v, v c ) (6)
with v c the closest feature point from v:

v c ∈ F/δ ′ (v, v c )=min v f i ∈F δ(v, v fi ) (7)
Notice that f m is invariant to uniform scaling (thanks to the normalization), rotation and translation (thanks to the use of geodesic distances). Figure 5 presents some computations of f m over arbitrary shapes, the number of extracted feature points (|F |)a ndt henumbe rofc r i t i c a lpoi nt s( |C|, identified according to the classification proposed in [START_REF] Cole-Mclaughlin | Loops in Reeb graphs of 2-manifolds[END_REF]). f m has been defined so as it tends to maxima (in green) at feature points and it tends to minima (in red) at the center of the object. As shown in figure 5, f m generates an important number of critical points. Consequently, standard Reeb graph construction algorithms would create large graphs, counting as many nodes as critical points, which is a major issue for topological skeleton extraction. In the next section, we present a formulation of discrete contours,whichenablesa unified graph construction and simplification process.

Discrete contours of a mapping function

In this section, we propose to construct discrete contours.I nt h en e x ts e c t i o n s ,t h o s ediscrete contours will be used either to detect topological changes or to detect curvature transitions, providing enhanced topological skeletons, without re-meshing and without any input parameter.

Defining contours of a real function f computed on a triangulated surface T is not a simple problem. In the continuous case, two points p 1 and p 2 belong to the same level line f -1 (f (p 1 )) if f (p 2 )f (p 1 )=0 .M o r e o v e r ,p 1 and p 2 belong to the same contour if they belong to the same connected component of f -1 (f (p 1 )).

In the discrete case, for a given vertex v ∈ T ,depending on T sampling, f -1 (f (v)) is often reduced to the vertex v itself. With regard to definition 1, a correct Reeb graph could not be constructed from this formulation of discrete contours, because the conditions of the equivalence relation would rarely be satisfied.

To preserve contour topological properties in the dis- crete case, we define the discrete level line Γ(v) associated to the vertex v as a curve computed along the edges of T which approximates by upper value the continuous level line f -1 (f (v)).F i g u r e6s h o w sdiscrete level lines traversing an arbitrary triangulation, with regard to the height function. Moreover, each connected subset of a discrete level line is referred as a discrete contour.I np a r t i c u l a r ,w ed efine the discrete contour γ(v) associated to the vertex v as the connected subset of Γ(v) containing v.N o t i c et h a tt h e more T will be dense, the more discrete contours will tend to continuous contours. Discrete contours can be computed for the whole mesh using a step by step gradient ascent process, described in algorithm 1. It handles two heaps, respectively the set of visited vertices Vtand the set of candidate vertices for visit Cd.A te a c hs t e p ,Cd surrounds Vtby upper value.

Algorithm 1 Discrete contour computation.

Vt=

∅ Cd ←{argmin v ′ ∈T f (v ′ )} while Cd = ∅ do v ← argmin v ′ ∈Cd f (v ′ ) 5: Γ(v) ← Cd γ(v) ← connected subset of Cd,containingv
Remove v from Cd Cd ← Cd ∪{v neighbors, which are not in Vt} Add v to Vt 10: end while A discrete level line locally minimizes its difference with the continuous level line it approximates. This difference can be expressed as follows:

vi∈Γ(v) (f (v i ) -f (v)) /f (v i ) ≥ f (v), ∀v i ∈ Γ(v) (8)
In expression 8, f (v) is a constant term. Consequently, minimizing this expression is equivalent to minimizing f (v i ) which is actually performed at each iteration of algorithm 1. Cd always surrounds Vt by upper value and minimizes expression 8, thus it is equivalent to Γ(v).

In figure 7, examples of discrete level lines Γ(v i ) are shown, at different iterations i of the algorithm. Vt vertex set is displayed in white and Γ(v) is displayed in red.

Visiting in a recursive fashion each vertex of Γ(v) enables the identification of each of its connected subsets, and particularly γ(v).

Topological analysis of discrete contours

Standard Reeb graph construction algorithms need simplification in order to remove noisy details. In this section, we propose a unified algorithm for graph construction and simplification, based on the topological analysis of discrete contours. Following the definition 1 of a Reeb graph in the continuous case, we can state an analog equivalence relation in the discrete case between two vertices v 1 ,v 2 ∈ T , based on our notion of discrete contour:

(v 1 ,f(v 1 )) ∼ (v 2 ,f(v 2 )) ⇐⇒ v 2 ∈ Γ(v 1 ) v 2 ∈ γ(v 1 ) (9)
v 1 and v 2 belong to the same connected component if they satisfy the above conditions. Therefore, at each iteration of the contour computation algorithm, each individual connected component of T ,traversedbyΓ(v),canbeidentified. Thus, topological changes can be detected observing the number N Γ(v) of connected subset of Γ(v),asf evolves. We define three types of topological changes:

1. bifurcations:w h e nN Γ(v) increases from iteration t to iteration t +1(Γ(v) splits in two contours from 8(a) to 8(b)), 2. junctions:w h e nN Γ(v) decreases from t to t +1and when several discrete contours merge (two contours merge in one from 8(b) to 8(c)), 3. terminations:w h e nN Γ(v) decreases from t to t +1, without discrete contour merge. Figure 9 shows several dual Reeb graphs obtained with this strategy, with regard to f m .C onne c t e dc ompone nt sa re represented by the nodes located at their barycenter.

The main contribution of our algorithm is that graph construction and simplification are performed at the same time. If we compare figures 9 and 5, we notice that the dual Reeb graphs do not reflect the presence of noisy critical points (points in red in figure 5), because discrete level lines do not disconnect in those configurations. Standard Reeb graph algorithms would have generated graphs counting as many nodes as critical points -94 nodes for the hand model and 92 nodes for the horse model -while in our approach only meaningful topological variations are encoded in the graph.

In comparison to [START_REF] Attene | Shape understanding by contour-driven retiling[END_REF][START_REF] Hilaga | Topology matching for fully automatic similarity estimation of 3D shapes[END_REF], no re-meshing and no input parameter, such as a slicing parameter, are required. Consequently, as no assumption is made about T sampling, this algorithm is robust against variations in mesh sampling, as shown in fig. 10. Furthermore, as f m is based on normalized geodesic distance evaluations, presented graphs are also invariant to geometrical transformations (rotation, translation and uniform scaling).

In this section, we presented a unified graph construction and simplification algorithm, based on the topological analysis of discrete contours.A sc o n t o u r sd on o td i s c o n n e c ti n f m noisy parts, resulting graphs reveal the shape most significant features. However, a strict topological analysis cannot discriminate visually meaningful sub-parts of a given connected component. To overcome this issue, we propose in the next section to analyze the geometrical characteristics of discrete contours to detect constrictions.

Geometrical analysis of discrete contours

Constriction approximations enable the subdivision of the branches of topological skeletons into more visually interesting parts. In this section, we propose an algorithm for constriction approximation, based on the analysis of the curvature of discrete contours.F o re a c hdiscrete contour identified in the previous stage, we compute its Gaussian curvature and we identify local minima as constrictions.

Model

f τ =3 f τ =5 f τ =10 f τ =15

Topological constraint

Since constrictions are defined as closed curves, the analysis has to be restricted on closed discrete contours only. Considering each contour γ(v) as a connected and nondirected planar graph G, γ(v) is a cycle, and consequently ac l o s e dc u r v e ,i ft h ed e g r e eo fa l li t sv e r t i c e se q u a l st w o . Therefore for each discrete contour of T reduced to a planar graph G,t h ed e g r e eo fe a c ho fi t sv e r t e xi sc o m p u t e d and we only consider in the rest of our algorithms contours that satisfy the above property.

Concavity curves

In our experiments, the average curvature ζ(γ(v)) of each discrete contour γ(v) is estimated by computing the Discrete Gaussian Curvature [START_REF] Meyer | Discrete differential-geometry operators for triangulated 2manifolds[END_REF] in each of its vertex. If ζ(γ(v)) is positive, γ(v) neighborhood is globally convex, otherwise it is concave. Constrictions appear on the narrowest, or the most concave, parts of a surface. Consequently, in order to only consider concave discrete contours, we compute ζ ′ (γ(v)) as follows:

ζ ′ (γ(v)) = ζ(γ(v)) if ζ(γ(v)) ≤ 0 0 if ζ(γ(v)) > 0 (10)
During the discrete contour computation, each contour is stored in the related node of the dual Reeb graph. As this algorithm visits T from the lowest to the highest f m values, for a given node of the graph, contours are automatically sorted by f m values. Computing ζ ′ (γ(v)) sequentially for each of these sorted contours gives, for a given node, a concavity curve,anovervie woftheconca vitye volutionasf m evolves.

Curves shown in figure 11 give examples of such evolutions, computed on the neck of the horse model (fig. 12(j)). The left values of these concavity curves correspond to the concavity estimations of the discrete contours located at the basis of the neck. The right values correspond to the concavity estimations of the discrete contours located at the end of the neck (basis of the head sub-part).

Constriction approximation

Curvature is a well-known noise sensitive entity. Consequently, to compute nice-looking approximations of constrictions, we have to reduce high frequency noise in concavity curves. Reducing noise on a one-dimensional data set is a trivial signal-processing problem. This can be achieved by applying an ideal low-pass filter of cutoff frequency f τ ,definedbythefollowingtransferfunction:

H(f γ(v) )= 1 if f γ(v) ≤ f τ 0 if f γ(v) >f τ (11)
Afi l t e r e dv e r s i o no fζ ′ (γ(v)) is given by the following expression, where FT stands for the Fourier Transform:

ζ ′ (γ(v)) = FT -1 (H(f γ(v) ) × FT(ζ ′ (γ(v)))) (12)
As shown in figure 11(b), low-pass filtering enables the discrimination of strongly concave contours. Consequently, for each node of the topological skeleton, we identify as constriction approximations the discrete contours that strongly minimize ζ ′ (γ(v)).T h e n ,t h ed u a lR e e bg r a p h is refined, subdividing each node using its constrictions as boundaries between sub-parts.

Experiments and results

In this section, we present and comment on experimental results obtained with our method and we discuss about its applications, particularly mesh deformation. Presented models are connected triangulated surfaces extracted from the Princeton Shape Benchmark database [START_REF] Shilane | The Princeton shape benchmark[END_REF].

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n) (o)

Discussion

Figure 12 presents intermediary results and enhanced topological skeletons of standard models.

Firstly, we can notice that our feature point extraction algorithm achieves correct extractions, even on complex areas, such as the hair of the humanoid model.

Secondly, our constriction approximation algorithm computes nice-looking constrictions (in red in fig. 12) even for coarsely designed objects (figs. 12(f) and 12(g)). It automatically adjusts its detection criterion to the connected component it is processing, enabling the identification of constrictions even on strongly tubular mesh sub-parts (like the legs of the horse models).

Table 1 shows that the number of identified constrictions is quite stable when f τ varies. The most visually interesting results have been obtained setting f τ =10.Thissettinghas been chosen for each model in fig. 12. However, constrictions are approximated assuming they appear along identified discrete contours only. This is a strong hypothesis. As shown in figure 12, thanks to our mapping function definition, this limitation is not detrimental when dealing with natural objects because some contours are actually identified on the articulations of prominent components.

Constriction approximation leads to the subdivision of the branches of the topological skeleton into more visually interesting parts: such as the decomposition of fingers into phalanxes (figs. 12(m) and 12(n)) or the decomposition of legs into thighs, calves and feet (figs. 12(k), 12(l) and 12(o)).

Thanks to f m invariance properties, those skeletons are invariant to geometrical transformations: translation, rotation and uniform scaling. Furthermore, as shown section 6, thanks to our graph construction strategy, no noisy details are encoded in the skeletons.

Time complexity

Given an input connected triangulated surface T ,l e t n be the number of vertices in T .F e a t u r e p o i n t e x t r a ction is performed in O 12. Notice that our overall method has a significantly lower running time than latest constriction detection [START_REF] Hétroy | Constriction computation using surface curvature[END_REF] or skeleton extraction [START_REF] Wu | Domain connected graph: the skeleton of a closed 3D shape for animation[END_REF] methods, for equivalently sampled meshes.

(n × log(n)) steps. f m is com- puted in O(|F |×n × log(n))

Example of application: mesh deformation

Topological skeletons have shown to benefit various applications in computer graphics [START_REF] Biasotti | An overview on properties and efficacy of topological skeletons in shape modelling[END_REF]. For example, within the framework of shape retrieval, thanks to their invariance properties, enhanced topological skeletons can be used for shape similarity estimation. They are also good supports for shape compression, metamorphosis, texture mapping, etc.

In this paper, to show the usability of our approach, we focus on mesh deformation. Each node of the enhanced topological skeletons references each vertex of the related mesh sub-component. Thus, a novice user can easily apply deformations on selected parts of the object.

Since mesh deformation is not in the scope of this paper, in our experiments, models are deformed by applying simple rotations to components, but more sophisticated strategies can be used [START_REF] Lewis | Pose space deformations: A unified approach to shape[END_REF]. Given an angle and an axis of rotation, a rotation matrix is computed. Then it is applied to each vertex of the selected node, providing nice-looking deformations, as shown in figures 1(f) and 13.

Conclusion and future works

In this paper, we presented a fully automatic algorithm for the extraction of affine-invariant enhanced topological skeletons. It first computes a dual Reeb graph. Then it refines it using constrictions as boundaries between mesh subparts. To our knowledge, this is the first approach which unifies Reeb graph and constriction computations.

Our scientific contribution resides in three points. Firstly, we proposed a robust and straightforward feature point extraction algorithm. It enables the computation of an invariant mapping function which reveals well the shape most significant features. Secondly, we presented an algorithm for discrete contours computation. We showed that a topological analysis of these discrete contours enables a unified Reeb graph construction and simplification process. Resulting graphs do not encode noisy details and they are robust against variations in mesh sampling and invariant to geometrical transformations. Finally, we showed that a geometrical analysis of the discrete contours provides nice-looking constriction approximations on prominent components, enabling the refinement of dual Reeb graphs into more visually meaningful skeletons.

Our algorithm computes skeletons with satisfactory execution times, without any input parameter or pre-processing stage. Consequently, it is a good candidate for various applications in computer graphics, like shape deformation (experimented in this paper), retrieval, metamorphosis, compression, texture mapping, etc.

In the future, we would like to experiment more robust geometrical analyses and constriction sliding algorithms,in order to provide more visually interesting skeletons. Moreover, forcing the position of the skeleton inside the object and preserving shape symmetry are enhancements which benefit certain applications and which are currently under investigation.

Figure 1 .

 1 Figure 1. Main steps of presented method on a standard model.

Figure 2 .

 2 Figure 2. Evolution of the level lines of the height function on a bi-torus, its critical points and its Reeb graph.

  (a) E 1 . (b) E 2 . (c) E 1 ∩ E 2 .

Figure 3 .

 3 Figure 3. Feature point extraction overview.

  (b)). Extremities of prominent components are configurations where f g1 and f g2 tend to an extremum (figs. 3(a) and 3(b)). Consequently, the set of feature points is both included in E 1 and E 2 .T h e r e f o r e ,w ed e fi n et h es e t of feature points F of T (fig. 3(c)) as follows:

Figure 4 .

 4 Figure 4. Feature point extraction robustness against mesh sampling variations.

  (a) |F | =2 , |C| =326. (b) |F | =6 , |C| =94. (c) |F | =7, |C| =92.

Figure 5 .

 5 Figure 5. f m mapping function computed on arbitrary shapes.

  (a) Γ(va). (b) Γ(v b ).

Figure 6 .

 6 Figure 6. Example of continuous (red) and discrete (green) level lines (height function).

Figure 7 .

 7 Figure 7. Example of discrete level lines on a 25 000 vertex mesh (f m function).

Figure 8 .

 8 Figure 8. Bifurcation and junction contexts on a torus shape (height function).

Figure 9 .

 9 Figure 9. Dual Reeb graphs of primitive and complex shapes (f m function).

Figure 10 .

 10 Figure 10. Algorithm robustness against mesh sampling variations.

  (a) Unfiltered curve. (b) Filtered curve (fτ =1 0 ).

Figure 11 .

 11 Figure 11. Concavity curves for the neck of the horse model (unfiltered and filtered).

  (a) Input mesh. (b) Enhanced topological skeleton. (c) Deformed skeleton. (d) Deformed mesh.

Figure 13 .

 13 Figure 13. Application to mesh deformation.

  steps with |F | the number of identified feature points. Notice that f m has a lower computational cost than the function proposed in [10] (|F | rarely exceeds 20). Each discrete contour computation takes O(log(n)+n).T h e r e f o r e ,a sc o n t o u r sa r ec o m p u t e df o r each vertex in T ,t h eo v e r a l ld i s c r e t ec o n t o u rc o m p u t ation takes O(n 2 ) steps. Topological and geometrical analyses are more straightforward. Topological analysis is performed in O(n) steps. Concavity curves are computed in O(n) steps. Their smoothing is realized in O(n × log(n)), using the Fast Fourier Transform algorithm. Consequently, we can state that the overall complexity of our method is bounded by the discrete contour computation, which takes O(n 2 ) steps. Presented algorithms have been implemented in C language under GNU/Linux and experimented on a desktop PC with a 3 GHz P4-CPU and 2 gigabytes of RAM. Table 2s h o w st h ec o m p u t a t i o nt i m e sc o r r e s p o n d i n gt ot h em o dels presented in fig.

Table 1 . Number of constrictions with differ- ent concavity curve cutoff frequencies (f τ ).

 1

	Humanoid	4	8	9	12
	Horse 1	11	11	11	12
	Hand 1	5	8	11	12
	Hand 2	6	9	11	11
	Horse 2	9	14	19	19

Figure 12. Feature points, constriction approximations andenhancedtopologicalskeletonsofstan- dard models.

	Model	Faces Feature pts. Constrict. Time
	Humanoid 1900	19	9	0.5 s.
	Horse 1	20 000	10	11	12 s.
	Hand 1	50 000	6	11	75 s.
	Hand 2	50 000	6	11	100 s.
	Horse 2	40 000	7	19	35 s.

Table 2 . Computation times.

 2


Julien Tierny
email: tierny@lifl.fr

Jean-Philippe Vandeborre

Mohamed Daoudi
email: daoudi@enic.fr

Invariant High Level Reeb Graphs of 3D Polygonal Meshes

a) Feature points. (b) Invariant quotient function. (c) Skeletal representation of the high level Reeb graph. (d) Application to mesh deformation.

Figure 1. Overview of presented method on an arbitrary mesh.

Introduction

Polygonal mesh is a widely used representation of 3D shapes. However, many applications in computer graphics need a higher level shape description as an input, like its structure for example.

To solve this problem, several approaches have been developed, like mesh segmentation [START_REF] Shamir | A formalization of boundary mesh segmentation[END_REF] or skeleton extraction [START_REF] Gagvani | Parameter controlled volume thinning[END_REF]. Topological methods, based on Reeb graphs, present the advantage to preserve the topological properties of the surface. Unfortunately, in practice, the construction of Reeb graphs raises several issues, like invariance constraint nonrespect or non-significant critical point identification. This may come to results of restricted usability and of low semantic interest, as underlined in [START_REF] Ni | Fair Morse functions for extracting the topological structure of a surface mesh[END_REF], that we refer as lowlevel Reeb graphs.

In this paper, we present a novel method for the construction of invariant high level Reeb graphs. First, we present theoretical backgrounds and related works. Secondly, we introduce a new feature point extraction algorithm (cf. figure 1(a)), which is used for the computation of a meaningful and invariant quotient function (figure 1(b)), and a new graph construction algorithm (figure 1(c)) that excludes non-significant critical points. Finally, we present and comment on experimental results and evoke possible applications, like mesh deformation (figure 1(d)).

Theoretical background

A Reeb graph [START_REF] Reeb | Sur les points singuliers d'une forme de pfaff complètement intégrable ou d'une fonction numérique[END_REF] is a structure that depicts the evolutions of the level lines of a given scalar function, usually a Morse function [START_REF] Morse | Relations between the critical points of a real function of n independant variables[END_REF], defined over objects of any dimension (k-manifolds). In this section, we introduce several notions, restricted to 2-manifolds, used in our approach.

Definition 1 (Level line) Let f be a scalar function defined on a 2-manifold M , f : M → R. Let L be the set of points p i , such as ∀p i ∈ L, f (p i ) = f L . L is the level line corresponding to the value f L , noted f -1 (f L ). Concretely speaking, a Reeb graph is composed of nodes, which represent f critical points (see illustrations 2 and 3), and of edges, which represent the connected components linking f critical points. Most of the Reeb graph construction algorithms first identify the set of f critical points and secondly build the graph from the connectivity relations of these points.

Definition 3 (Critical point) Let f be a scalar function defined on a

2-manifold M , f : M → R. A point p ∈ M is a critical point of f is the gradient of f vanished in p.
Definition 4 (Non-degenerate critical point) A critical point p of a scalar function f defined on a 2-manifold is called a non-degenerate critical point if the matrix of second partials is non-singular in p. In the discrete case, given a vertex v of a triangulated surface T , let L + and L -be respectively the number of connected sets of vertices v i and v j on Lk(v) (v direct neighbors, its link) such as f

(v i ) > f (v) and f (v j ) < f (v).
Vertices corresponding to non-degenerate critical points can be classified into three categories, according to L + and L -, as shown in figure 3: local minima (3(a)), local maxima (3(b)) and simple saddles (3(c)). We will use figure 3 color convention in the rest of the paper.

As mentioned by Biasotti et al. [START_REF] Biasotti | An overview on properties and efficacy of topological skeletons in shape modelling[END_REF], continuous functions, and particularly Morse functions [START_REF] Morse | Relations between the critical points of a real function of n independant variables[END_REF] (when all the critical points are non-degenerate), are pertinent choices for functions f . For clarity puprose, we will refer to candidate functions as quotient functions in the rest of the paper, to underline their use in graph construction.

Figure 2 illustrates previous definitions, presenting respectively a scalar function (the height function) computed on a 2-manifold, its non-degenerate critical points, level lines of arbitrary range and the corresponding Reeb graph.

Related work

Several approaches have been explored by the computer graphics community in polygonal mesh decomposition.

Topological methods aim to produce structural descriptions of shapes where identified sub-components form connected sets. Morse and Reeb graph theories are two powerful mathematical tools that respectively enable to identify topological points of interest over the mesh and to capture their connectivity relations into a graph structure (see figure 2). Several algorithms [START_REF] Tarasov | Construction of contour trees in 3D in O((n) log (n)) steps[END_REF][START_REF] Carr | Computing contour trees in all dimensions[END_REF] have been proposed to construct Reeb-like graphs in O(n × log(n)) steps (where n is the number of vertices) from piecewise linear functions computed over the mesh. Tarasov and Vyalyi [START_REF] Tarasov | Construction of contour trees in 3D in O((n) log (n)) steps[END_REF] propose local mesh redefinition in order to purge multiple saddles and to benefit from Morse functions properties [START_REF] Morse | Relations between the critical points of a real function of n independant variables[END_REF]. However, those papers do not address the problem of the quotient function definition. Lazarus and Verroust [START_REF] Lazarus | Level set diagrams of polyhedral objects[END_REF] introduces a novel notion of topological structure called Level Set Di-agram, whose construction is somewhat similar to that of a Reeb graph. An interesting point of this method is the use of a quotient function, exploited for the diagram construction, based on the geodesic distance to a source vertex. Unfortunately, this scalar function suffers from a high instability [START_REF] Hilaga | Topology matching for fully automatic similarity estimation of 3D shapes[END_REF], due to the non-determinist election of the source vertex, which excludes its use in applications where stability is fundamental, like shape retrieval for example. To overcome this problem, within the framework of shape retrieval, Hilaga et al. [START_REF] Hilaga | Topology matching for fully automatic similarity estimation of 3D shapes[END_REF] propose to integrate this function over the whole mesh. More precisely, they define an approximation function that introduces geodesic distances to a set of base vertices, homogeneously spread over the mesh. Moreover, they propose a new topological structure scheme called Multiresolutional Reeb Graph. Basically, they proceed in a succession of mesh partitionings, according to different value-ranges of their quotient function. This initiative is of a major interest for 3D shape retrieval but does not fit other applicative contexts mainly because it will be difficult to automatically define a proper value-range for an application or another. For an interesting survey on computational topology approaches for shape modeling, we defer the reader to [START_REF] Biasotti | An overview on properties and efficacy of topological skeletons in shape modelling[END_REF].

In short, the first key issue of topological approaches resides in the definition of the quotient function. A non pertinent function will present an important number of critical points, as underlined in [START_REF] Ni | Fair Morse functions for extracting the topological structure of a surface mesh[END_REF]. Consequently the resulting Reeb graph will present many nodes and edges and will not afford a global and meaningful description of the shape (low-level Reeb graphs). Moreover, depending on application needs, the quotient function should present stability and invariance properties [START_REF] Lazarus | Level set diagrams of polyhedral objects[END_REF][START_REF] Hilaga | Topology matching for fully automatic similarity estimation of 3D shapes[END_REF] but those functions often generate many critical points from our experience.

The second key issue of topological approaches resides in the graph construction strategy, which should only use significant critical points and should not be conditioned by a user-parameter. In this paper, we present a novel method, whose objective is to produce high-level, stable, invariant and generically-exploitable Reeb graphs, solving previously mentioned issues.

Invariant high level Reeb graphs

Reeb graphs give a good overview of the structure of polygonal meshes. Unfortunately, constructing a Reeb graph from a non pertinent quotient function often leads to results of low semantic interestlow-level Reeb graphsbecause of the number of identified critical points.

In this paper, we propose to challenge this key issue by proposing an invariant quotient function combined with a pertinent critical point selection algorithm, that enables the construction of high-level Reeb graphs over piecewise continuous functions.

Invariant quotient function

Several scalar functions have been proposed by the computer graphics community to construct topological structures. The choice of this function will directly condition the stability properties of the topological structure and therefore its usability in application fields. In this section, we detail our strategy to compute an invariant and meaningful quotient function: first, we extract mesh feature points and secondly, for each vertex in the mesh, we compute its geodesic distance to the closest feature point.

Metric space

Before defining the scalar function itself, a metric space needs to be defined. We will distinguish two kinds of metric spaces: on the first hand euclidean ones and on the other hand Riemannian ones.

Rotation and translation invariance properties can be obtained in both, as long as origins are taken relatively to the mesh. As an example, with an euclidean metric space, the euclidean barycenter of the mesh can be chosen.

Scaling invariance property can be obtained in both metric spaces as well, normalizing values in a regard to mesh global extrema.

In shape modeling, geodesic distances are interesting Riemannian metrics, because their evaluations are tolerant to variations in model pose. In metric geometry, a geodesic is a curve which is everywhere locally a distance minimizer. More concretely, the geodesic distance between two vertices is the length of the shortest path between them along the mesh, according to a given metric. As an example, a geodesic distance from a hand of a humanoid model to its head will be the same whether its arms are folded or not. For a formal description of Riemannian geometry, we defer the reader to [START_REF] Fomenko | Topological Modeling for Visualization[END_REF].

In our approach, to guarantee invariance to rotation, translation, scaling and model pose, we decide to choose Riemannian metrics, and particularly geodesic distances. From an algorithmic point of view, geodesic distances can be approximated by the Moore-Dijkstra algorithm (distance minimizing in weighted graphs). In the rest of this paper, we will refer to δ(v 1 , v 2 ) as the normalized approximation of the geodesic distance from vertex v 1 to vertex v 2 .

Feature point extraction

Feature points are mesh vertices located on extremities of prominent components. From a perceptive point of view, the set of feature points describes in a meaningful way the global structure of a shape. From this observation, in a comparable way to [START_REF] Brostow | Novel skeletal representation for articulated creatures[END_REF], we propose in our method to use the set of feature points as origins for geodesic distance evaluations.

Several algorithms have been proposed for feature point extraction. They find applications in various fields, like shape metamorphosis, deformation transfer, texture mapping, etc. For example, Mortara and Pantanè [START_REF] Mortara | Affine-invariant skeleton of 3D shapes[END_REF] propose to select as features points the vertices where gaussian curvature exceeds a given threshold. Unfortunately, this method can miss feature points because of the threshold parameter and cannot resolve extraction on constant curvature areas. Katz et al. [START_REF] Katz | Mesh segmentation using feature point and core extraction[END_REF] developed an algorithm based on multidimensional scaling, in quadratic execution complexity.

Here, we propose a quite straightforward algorithm, based on topological tools. Most of the geodesic based scalar function local extrema appear at extremities of prominent components (see illustrations 4(a) and 4(b)), mainly because gradients vanish in those configurations (cf. definition 3). Therefore, we propose to realize a crossed analysis, using two geodesic based scalar functions -whose origins are the mesh most distant vertices -and to intersect the sets of their local extrema.

Definition 5 (Feature points) Let f 1 and f 2 be two scalar functions defined on a connected triangulated surface T as follows:

f 1 (v) = δ(v, v s1 ) (1)
with:

v s1 ∈ T / δ(v s1 , v r ) = max v∈T δ(v, v r ) (2)
with v r ∈ T a randomly chosen vertex, and:

f 2 (v) = δ(v, v s2 ) (3)
with:

v s2 ∈ T / f 1 (v s2 ) = max v∈T f 1 (v) (4)
Let E 1 and E 2 be the sets of local extrema of f 1 and f 2 . We define the set of feature points F of T as the intersection of E 1 and E 2 :

F = E 1 ∩ E 2 (5)
Concretely, we perform a crossed analysis in order to purge non-isolated extrema, as illustrated in figure 4: f 1 local extrema are displayed in yellow (figure 4(a)) and f 2 extrema in cyan (figure 4(b)). Equation 5 is well illustrated in figure 4.

In practice, f 1 and f 2 isolated local extrema do not appear exactly on the same vertices but in the same geodesic neighborhood. Therefore, we relax the intersection constraint as follows: In our experiments, setting ǫ = 0.05 gives accurate results. With this configuration, m never exceeds two percent of n and k rarely exceeds 30 (depending on the model's topological complexity). Therefore E 1 and E 2 relaxed intersection algorithm's execution complexity is negligible compared to f 1 and f 2 computations.

v ∈ F ⇐⇒        ∃v e1 ∈ E 1 / δ(v, v e1 ) < ǫ ∃v e2 ∈ E 2 / δ(v, v e2 ) < ǫ δ(v, v fi ) > ǫ ∀v fi ∈ F ǫ ∈ [0, 1] (6) (a) E 1 . (b) E 2 . (c) E 1 ∩ E 2 .
In this section, we presented a straightforward algorithm for mesh feature point extraction. This algorithm is based on geodesic distance evaluations and therefore is stable and invariant to geometrical transformations and to model pose. Moreover, in order to select feature points, we observe geodesic gradient behaviors. Consequently, we can state that our method is robust against mesh sampling variations, as illustrated in figure 5.

Quotient function definition

The quotient function definition depends on what is expected to be revealed. As an example, for terrain modeling, the height function will present critical points over hills and valleys and will afford consequently an appropriate topological description. In our approach, we would like to identify meaningful mesh sub-components.

As underlined in paragraph 4.1.2, the set of feature points describes in a meaningful way the global structure of a shape. This set is invariant to geometrical transformations and variations in model pose and mesh sampling. From these observations, we propose to compute a meaningful quotient function from geodesic distances to the set of feature points.

In particular, an interesting objective would be to make f function level lines cut as precisely as possible the basis of prominent components, to afford a meaningful decomposition. With this aim, we use the set of feature points as origins for distance evaluations and propose the following quotient function, computed for each vertex v of an input connected triangulation T , noted f q in the rest of the paper:

f q (v) = 1 -δ(v, v c ) (7)
with v c the closest feature point from v:

v c ∈ F / δ(v, v c ) = min v f i ∈F δ(v, v fi ) (8)
Figure 6 presents some computations of f q over arbitrary shapes, as well as the number of extracted feature points.

In our experiments, we compute the set of geodesic distances towards feature points within the feature point extraction algorithm. Therefore, the number of iterations of the Moore-Dijkstra algorithm over the whole mesh is equal to the number of feature points. The execution complexity of f q computation can be consequently approximated by O(|F | × n × log(n)), with n the number of vertices in the mesh.

High level Reeb graph construction

As we can see in figure 6, f q is a piecewise continuous function only. More precisely, discontinuity appears on areas where geodesic origins change. In figure 6, those areas can be identified with the rings of local minima and local saddles (in red and black). From a theoretical point of view, resulting Reeb graph should present as many nodes as function critical points and consequently would not afford a global, or high level, description of the shape.

In this section, we propose to challenge this issue with a novel critical point election algorithm. With this aim, we propose to focus on an intuitive description of Reeb graphs: Reeb graphs are topological structures that depict the connectivity evolution of the level lines of a given scalar function defined on a 2-manifold.

Consequently, we present a method that does not propose to construct a Reeb graph from the connectivity relations of f critical points. Instead, we propose to construct discrete level lines around f global minimum and to observe their connectivity relations while f evolves.

Discrete level lines evolution can be modeled as a geodesic propagation within a metric space based on the quotient function values. Such a propagation can be obtained with the Moore-Dijkstra algorithm, considering for each vertex its f value as weight.

Let t be a parametric variable that denotes an iteration of the geodesic propagation algorithm. Now, we can introduce the notion of discrete level line:

Definition 6 (Discrete level line) Let f be a scalar function defined on a connected triangulated surface T , f :

T → R.
Let C(t) be the set of candidate vertices for absorption at the iteration t of the geodesic propagation algorithm.

Let l 1 (t), l 2 (t), ..., l k (t) be the connected subsets of vertices belonging to C(t). We define each connected subset l i (t) as a discrete level line. Figure 7 gives an overview of the geodesic propagation in the metric space defined by f q . The area depicted in white denotes the vertices visited by the Moore-Dijkstra algorithm, whereas the red areas correspond to the vertices belonging to C(t). Moreover, each connected subset of C(t) (red thin lines) is referred as a discrete level line.

This notion of discrete level line can be compared to the one of topological ring, presented by Mortara and Pantanè [START_REF] Mortara | Affine-invariant skeleton of 3D shapes[END_REF]. However, their graph construction strategy is highly dependent on their quotient function definition. Moreover, to benefit from the properties of geodesic distances, their quotient function computation needs a regular mesh refinement, which increases the overall time complexity of the approach and decreases its tolerance to mesh sampling variations.

In our algorithms, we analyze the connectivity evolutions of the discrete level lines at each iteration t of the geodesic propagation. Moreover, we construct simultaneously a dual Reeb graph (where connected components are represented with nodes and connectivity relations with edges), according to the following notions of topological variations: bifurcations, junctions and terminations. Definition 7 (Bifurcation) Let f be a scalar function defined on a connected triangulated surface T . Let L(t) = {l 1 (t), l 2 (t), ..., l k (t)} be the set of discrete level lines l i (t) at iteration t of f propagation and v t ∈ T the last visited vertex. A bifurcation happens in v t iff:

|L(t)| > |L(t -1)| (9)
Figure 8 presents the bifurcation and junction contexts on a torus with the height function. In 8(a), L(t) is composed of only one discrete level line, which splits in two in 8(b): a bifurcation is created on the graph (figure 8(e)). Similarly, we introduce the notion of junction as follows: Definition 8 (Junction) Let f be a scalar function defined on a connected triangulated surface T . Let L(t) = {l 1 (t), l 2 (t), ..., l k (t)} be the set of discrete level lines l i (t) at iteration t of f propagation and v t ∈ T the last visited vertex. A junction happens in v t iff:

|L(t)| < |L(t -1)| ∃v n ∈ Lk(v t ) / v n ∈ C(t) (10)
In 8(b), L(t) is composed of two discrete level lines, which merge in one in 8(c): a junction is created on the graph (figure 8(f)). Finally, we introduce the notion of termination: Definition 9 (Termination) Let f be a scalar function defined on a connected triangulated surface T . Let L(t) = {l 1 (t), l 2 (t), ..., l k (t)} be the set of discrete level lines l i (t) at iteration t of f propagation and v t ∈ T the last visited vertex. A termination happens in v t iff:

|L(t)| < |L(t -1)| v n / ∈ C(t), ∀v n ∈ Lk(v t ) (11)
At each step t of the geodesic propagation, we reconstruct L(t) and apply needed topological variations on the graph (figures 8(e) and 8(f)), according to |L(t)| evolutions (equations 9, 10 and 11), as illustrated in the next section.

Experimental results and comments

In this section, we present and comment on experimental results obtained with our method and discuss about its applications. Presented models are connected triangulated surfaces extracted from the Princeton Shape Benchmark database [START_REF] Shilane | The Princeton shape benchmark[END_REF].

Time complexity

Let n be the number of vertices in the input mesh. As mentioned section 4.1.2, the feature point extraction is realized in O(n × log(n)) steps. Moreover, f q is computed within the feature point extraction. Consequently, the invariant quotient function computation time complexity is bounded by the Moore-Dijkstra algorithm and is realized in

O(|F | × n × log(n)).
As for the graph construction algorithm, we model C(t) with a binary priority heap, which means that addition and deletion of vertices are performed in O(log(n)). At a given iteration t of the geodesic propagation, in order to observe topological variations defined in 4.2, we re-construct each discrete level line with a recursive algorithm, in O(n) steps. Consequently, an iteration t of the geodesic propagation is realized in O(n + log(n)) steps. Therefore, as the complete geodesic propagation takes n iterations, the overall complexity of the graph construction takes O(n 2 ) steps.

Immersion strategies

By definition, a Reeb graph is a graph representation only. Therefore, for many applications, it is mandatory to define an R 3 immersion strategy.

In this paper, for illustration purpose, we use the following strategy: each node of the high level Reeb graph, corresponding to a whole connected component, is placed at the euclidean barycenter of its related component, as shown in figure 9. Notice that the root component, with bigger radius, represents the neighborhood of f global minimum. This area is isolated when the geodesic wavefront first collapses on itself, at the beginning of the algorithm.

Topological skeletons are a Reeb graph variant that can be obtained by placing a point at the euclidean barycenter of each discrete level line, as shown in figure 10(a). Such a skeleton is a particularly pertinent shape descriptor for automatic mesh animation. With this aim, since a point of the topological skeleton references a discrete level line, it is possible to move each vertex of the related level line with a parameterized displacement, so as to get a smooth movement of the whole connected component. In our experiments, we deformed models by applying recursively simple rotations to components, as shown in figure 10(b), but more sophisticated strategies can be used, like in [START_REF] Katz | Hierarchical mesh decomposition using fuzzy clustering and cuts[END_REF].

High level Reeb graph properties

Stability of the high level description First, as illustrated in figure 9, both complex and primitive shapes are well managed by our algorithm. Moreover, it can handle non-null genus surfaces correctly: in figure 9(b), the hair of the humanoid model merges with its back which forms consequently a cycle in the graph. Secondly, we can say that the high level Reeb graphs of models belonging to the same class (figure 9(b), 9(c), 9(d) and 9(e)) are quite similar, which denotes the stability of our algorithm.

Finally, our geodesic propagation within the metric space defined by the quotient function, combined with our critical point election strategy, leads to accurate results: none of the presented graphs reflects the presence of critical points on discontinuity areas. With traditional Reeb graph construction algorithms, the graph of the model 9(f) would have counted about ninety critical points (see critical points in figure 6(b) for comparison). With our approach, only meaningful topological variations are encoded.

Affine invariance Thanks to our choice of metric space (geodesic distances, with origins taken relatively to the mesh), it is quite obvious that our approach is invariant against geometrical transformations: translation, rotation and uniform scaling. Robustness to mesh sampling Thanks to our notion of discrete level line, no hypothesis about mesh sampling is required. Therefore, we can state that our Reeb graph construction algorithm is tolerant against variations in mesh sampling, as shown in figure 11.

Thanks to those properties, high level Reeb graphs can be used in applications where invariance is fundamental. In shape retrieval for example, those graphs can be extended with geometrical attributes and shape comparison can be achieved with graph matching algorithms [START_REF] Hilaga | Topology matching for fully automatic similarity estimation of 3D shapes[END_REF].

Conclusion and future works

In this paper, we introduced a novel method for the construction of invariant high level Reeb graphs, topological entities that afford a global understanding of shapes, with satisfactory execution times and without input parameters. This method is composed of three main steps. First, we extract feature points thanks to a robust and straightforward algorithm. Then, we use feature points as geodesic origins for the computation of an invariant quotient function f q , used for topological analysis. Finally, we developed a new graph construction algorithm, based on a geodesic propagation in the metric space defined by f q . It observes the connectivity evolutions of discrete level lines and provides meaningful Reeb graphs, which only encode significant topological variations. We illustrated the utility and the accuracy of our approach with an application to mesh deformation.

The invariance properties of presented graphs (geometrical transformations, model pose and mesh sampling) make them good candidates for various applications in computer graphics, like shape animation, retrieval, compression, etc.

In the future, we would like to refine our mesh decomposition scheme in order to propose high level shape descriptions of higher semantic interest.

Figure 2 .

 2 Figure 2. Evolution of the level lines of the height function on a bi-torus, its critical points and its Reeb graph.

Definition 2 (

 2 Reeb graph) Let f : M → R be a scalar function defined on a 2-manifold M . The Reeb Graph of f is the quotient space of f in M × R by the equivalence relation (p 1 , f (p 1 )) ∼ (p 2 , f (p 2 )), verified if and only if:    f (p 1 ) = f (p 2 ) p 1 and p 2 belong to the same connected component of f -1 (f (p 1 ))

Figure 3 .

 3 Figure 3. Non-degenerate critical point classification.

Figure 4 .

 4 Figure 4. Feature point extraction overview.

Figure 5 .

 5 Figure 5. Feature point extraction robustness against mesh sampling variations.

  50 000 vertices, 7 features points.

Figure 6 .

 6 Figure 6. Evolution of the level lines of f q and its critical points on arbitrary shapes.

Figure 7 .

 7 Figure 7. Geodesic propagation overview on a 25 000 vertex mesh (f q function).

Figure 8 .

 8 Figure 8. Bifurcation and junction contexts on a torus shape (height function).

Figure 9 .

 9 Figure 9. High level Reeb graphs of primitive and complex shapes.

Figure 10 .

 10 Figure 10. Skeletal representation of a high level Reeb graph (a) and an application to mesh deformation (b).

Figure 11 .

 11 Figure 11. High level Reeb graph construction algorithm robustness against mesh sampling variations. Robustness to model pose Robustness to model pose can be observed in figures 9(c) and 9(d): the graphs are similar wether the arms of the humanoid are folded or not.Robustness to mesh sampling Thanks to our notion of discrete level line, no hypothesis about mesh sampling is required. Therefore, we can state that our Reeb graph construction algorithm is tolerant against variations in mesh sampling, as shown in figure11.Thanks to those properties, high level Reeb graphs can be used in applications where invariance is fundamental. In shape retrieval for example, those graphs can be extended with geometrical attributes and shape comparison can be achieved with graph matching algorithms[START_REF] Hilaga | Topology matching for fully automatic similarity estimation of 3D shapes[END_REF].