Marcin Derlatka
email: m.derlatka@pb.edu.pl

Modified kNN algorithm for improved recognition accuracy of biometrics system based on gait

Keywords: k-nearest neighbors, human gait recognition, ground reaction forces

The k-nearest neighbors classifier is one of the most frequently used. It has several interesting properties, though it cannot be used in utilitarian biometric systems. This paper proposes the modification of kNN algorithm which ensures correct work even in the case of an attempt at access to the system by unauthorized people. Work of the algorithm was tested on data presenting ground reaction forces (GRF), generated during human's walk, obtained from measurements carried out on over 140 people. Differences between particular strides were determined with the Dynamic Time Warping (DTW). The recognition precision obtained was as high as 98% of biometric samples.

Introduction

Biometrics understood as a science of recognizing a human on the basis of his/her physical and/or behavioral features becomes increasingly popular and increasingly courageously enters our life. The number of devices (e.g. computers) and institutions (e.g. banks) where a particular type of biometrics is used is still increasing. Every biometric system, irrespective of the feature this system is based on, includes attribution of a 'biometric signature' provided by a person applying for access to resources, to one of the classes registered in the data base. Obviously, if the provided biometric signature is not similar enough the class pattern stored in the data base, there should appear information on a denial of access to the resources.

Among the most intuitive classifiers are the classifiers in which similarity is defined as the inverse of the distance in the feature space, e.g. the kNN classifier. In this case the input pattern is attributed to the class which includes most of its nearest neighbors. Intuitiveness, quite good generalizing properties and easiness of implementation result in the fact that the biometric systems working the basis of the kNN algorithm are frequent in literature [START_REF] Orr | The Smart Floor: a Mechanism for Natural User Identification and Tracking[END_REF], [START_REF] Tappert Ch | Performance of a Long-Text-Input Keystroke Biometric Authentication System Using an Improved K-Nearest-Neighbor Classification Method[END_REF]. Unfortunately, this classifier, despite its several modifications enhancing both the properties and the quality of classification [START_REF] Garcia -Pedrajas | Boosting K-nearest Classifier by Means of Input Space Projection[END_REF] [START_REF] Tappert Ch | Performance of a Long-Text-Input Keystroke Biometric Authentication System Using an Improved K-Nearest-Neighbor Classification Method[END_REF], has some features that prevent its direct application in biometric systems. They are:

attributing a given point of the space of features to one of the classes even when the distance, though minimal, is as long as we should talk about the lack of similarity to all patterns; in this situation the desired answer of the classifier is the information of unrecognizing the given person; while making decisions taking into consideration the constant number k of the nearest neighbors, even if part of them are relatively far from the point under classification.

These features, occurring also among other classifiers, lead to the fact that papers connected with biometrics very often assume that we deal with a closed group of people and we are not endangered with an attempt at getting to the resources by the people from outside of this group [START_REF] Yun | User Identification Using Gait Patterns on UbiFloorII[END_REF]. Certainly, such an approach limits the potential application field of a given system. One of more interesting behavioral biometrics is human gait. Gait is a very complex human activity. It is a symmetrical and repetitive phenomenon in its normal form. The movement of each limb is partitioned for the support phasewhen the foot is in contact with the ground and the swing phase when the foot is lifted and moved forward. The human gait is developed during child growth. It is assumed that the human gait pattern is evaluated till the child is seven years old and next stay almost unchanged till his death. The human gait is the result of synergistic activity of: bones, muscles and nervous systems. The cooperation of those three systems makes the human movement unique for every person.

Human recognizing using gait is not a new issue. It was Cutting and Kozlowski [START_REF] Cutting | Recognizing Friends by Their Walk: Gait Perception Without familiary Cues[END_REF] who already conducted researches which demonstrated that a human has a skill of recognizing the people he knows from a long distance by their way of moving, even if they wear clothes different than usual and have changed their hairstyle. Gafurov [START_REF] Gafurov | Used Authentication Based on Foot Motion[END_REF] distinguish the recent methods in human gait recognizing depending on the signals registered. They are methods based on the data obtained from:

video; floor sensor; portable sensors.

In case of video based methods the picture registered is usually converted frame by frame into silhouette sequences. Subsequently, depending on the applied methods there occurs an attempt at reading of selected parameter's of the human's gait and classification of a person [START_REF] Balista | Compact Time-independent Pattern Representation of Entire Human Gait Cycle for Tracking of Gait Irregularities[END_REF], [START_REF] Katiyar | A Study on Existing Gait Biometrics Approaches and Challenges[END_REF]. Advantages of these methods are undoubtedly a possibility of free motion of the person under examination, identification of many people simultaneously and these people's approval of the presence of cameras in buildings. Problems with recognizing people in the systems based on video cameras are generated by the examined person's change of clothes, a possibility of covering the examined person by objects or other people, changes in lighting as well as the sensitivity of certain parameters sought for to the position angle of the object in relation to the camera [START_REF] Katiyar | A Study on Existing Gait Biometrics Approaches and Challenges[END_REF]. These disadvantages do not occur in the group of methods based on the measurement of the interaction of the ground with the lower limbs of the person under examination (floor sensor based). Here the examined person must walk along an appropriately-prepared measurement path equipped with force plates [START_REF] Derlatka | Human Gait Recognition Based on Signals from Two Force Plates[END_REF] or a special floor with a network of photo interrupter sensors [START_REF] Yun | User Identification Using Gait Patterns on UbiFloorII[END_REF]. In the last method mentioned above the person under examination fully cooperates during the time of measurement. He/she is provided with measurement equipment, such as: accelerators [START_REF] Gafurov | Used Authentication Based on Foot Motion[END_REF], or opto-reflective markers [START_REF] Lin | Human recognition based on kinematics and kinetics of gait[END_REF].

The main objective of this paper is to demonstrate a modified kNN algorithm which can be applied in biometric systems. The work of the algorithm has been tested in author own research on human gait.

Ground Reaction Forces

In the biomechanical approach, the ground reaction force (GRF) is the force which is acting on the body as a response to its weight and inertia during the contact of the human plantar with the surface. The all three components of GRF are used in the presented work. They are: anterior/posterior F x , vertical F y and medial/lateral F z components of GRF. The common profiles of the GRF components are presented in Fig. 1 (a-c). The anterior/posterior component has two main phases. The value of F x is negative in the first phase. It is a result of the deceleration of the investigated lower limb, in this case the force direction is opposite in direction of walking. The minimum of the deceleration phase is most often reached a moment before the maximum of the limb-loading phase in the vertical component of GRF. The value of F x is positive in the second phase, respectively. The maximum of the acceleration phase is reached when the toe-off phase starts. There are three extremes in Fig. 1b. They correspond to: the maximum of the limb-loading phase, the minimum of the limb-unloading phase and the maximum of the propulsion phase (a moment before the toe off). It is not difficult to point to the same extremes for the medial/lateral component of GRF as for the vertical GRF.

The Modified k-NN Classification Algorithm

The k-nearest neighbors algorithm was modified for the needs of biometrics in such a way that it enables to make a decision depending on the degree of similarity between the case considered and the prototype patterns connected with a concrete user. To determine a similarity between the patterns author used a well-known DTW (Dynamic Time Warping) transform, which reduces the difference between time courses to the cost of adjusting one time series to another. The distance between n and m patterns has been calculated according to following formula:

D(n, m) = 6 c=1 D c (1)
where: D c is the DTW distance between the c components of GRF for patterns n and m. Values of similarity thresholds were defined for each user separately taking into consideration different dispersion of prototypes in the feature space for particular users. Consequently, for i-th user an average distance ρ i between particular prototypes was calculated. The distance limit connected with a concrete threshold was calculated according to formula:

ρ θi = (1.5 -θ)ρ i (2)
where: θ threshold; θ={0.1, 0.2, . . . , 0.9} and 1.5 is constant chosen arbitrary. It is worth noting that adopting such a definition means that for the threshold θ = 0.5 all the points whose distance from the prototypes is longer than the value of average ρ i of a given person will be treated as too dissimilar to a given class and will not be taken into account while making decisions (see point 4 of the following algorithm). The modified k-NN classification algorithm was as follows:

1. Determine distance D of the biometric pattern under examination from all prototype patterns in the data base. 2. Select k prototype patterns whose distances D to the pattern under examination are the shortest.

3. By majority vote determine the ID of the user in the data base. If two or more users are equally numerously represented among the patterns selected in point 2 (or the remaining after rejection in point 4) -select the one whose average distance from the selected pattern is the shortest. 4. Reject k' prototypes, for which distance D is longer than ρ θi for the given threshold θ of the user whose ID was selected in point 3. 5. Check in compliance with the procedure in point 3, whether ID for K=k-k' of prototypes remained unchanged. If so, we finish the classification assigning the examined biometric pattern to the ID class. If not, return to point 3. 6. In the case when k=k' (K=0) we recognize that a given biometric pattern cannot be classified in any of the classes at the assumed threshold θ.

Such a construction of the kNN classifier enables proper classification even if near the considered i-th class point in the input space there are more prototypes appertaining to competitive (j-th) class. This is possible on condition that they are more distant from the biometric signature under consideration than the limit threshold value for j-th class and the adopted threshold θ. During the examinations the person walked at a free speed in his/her own sport shoes on the measurement path, in which 2 Kistler's force plates were hidden, working at the frequency of 1 kHz. The volunteers performed several walks (14-20), as a result of which over 2500 strides were registered. 1056 patterns obtained from 132 people (57 man and 75 women) made a learning set (prototype points). 1500 patterns obtained from all the people under examination were treated as a testing set. The people under examination were divided into two groups: the so-called users or genuine of the system (132 people), who, as assumed, were to receive access to the system and the so-called impostors (10 people, 176 strides), who were represented only in the testing set and use to test the resistance of the biometric system to attempts at cheating (attempts at access to the data by unauthorized people).

For the obtained time courses we determined distances between all prototypes and all patterns tested in accordance with formula [START_REF] Balista | Compact Time-independent Pattern Representation of Entire Human Gait Cycle for Tracking of Gait Irregularities[END_REF]. In this paper we conducted a classification in accordance with the modified kNN algorithm proposed before, whereas for comparison we employed the classical version of the k-nearest neighbors. In both cases we assumed k=5. In the case of the classical kNN algorithm, the situation where among k-nearest neighbors was equally numerous representation of 2 or more classes, was treated as FRR error. The distances between the prototypes and the courses appertaining to the testing set performed by the same person were usually lower by a rank than in the case of the distances to these prototype points of strides of other people. The analysis of the situations in which an FRR error occurred in the classical kNN demonstrated that in the case of 19 strides it is possible to obtain a correct classification provided the most popular method of solving 'draws' were applied. It is enough, while making decisions, to adopt a criterion of average distance of the point under consideration to the prototypes of particular classes. Obviously, the remaining 10 strides will be classified incorrectly, which will raise the value of the FAR error up to 4.76%. It is important to note that the value of FAR error for the classical kNN is relatively low. Especially if we compare it to the results of other authors' works basing on signals of the same or similar types. Thus, for example, in the paper [START_REF] Orr | The Smart Floor: a Mechanism for Natural User Identification and Tracking[END_REF] (selected features of GRF and kNN profiles) the error amounted 7% on the group of 15 people, and in [START_REF] Yun | User Identification Using Gait Patterns on UbiFloorII[END_REF] (UbiFloorII and MLP): 1% and 10 people respectively. Slightly better results than in this paper were obtained in [START_REF] Pataky | Gait Recognition: Highly Unique Dynamic Plantar Pressure Patterns Among 104 Individuals[END_REF], where for foot pressure patterns, in the best cases, the error was 0.6% to 7.1% (104 people; 520 strides) and in [START_REF] Moustakidis | Feature Extraction Based on a Fuzzy Complementary Criterion for Gait Recognition Using GRF Signals[END_REF], where for GRF and so-called wavelet packet decomposition scheme the values of the error obtained were started from 0.5% (40 people).

It is important to underscore that this paper is based on much richer research material, which obviously affects the obtained results. In the case of the proposed modified kNN algorithm the value of FAR errors for the users independently from the adopted threshold are considerably lower than for the classical version of kNN. The analysis of the situations where the classical kNN generated an error allowed for demonstrating in which cases the proposed kNN algorithm works better. One of such cases occurs when a correct class has less numerous representation than a competitive class among the selected k prototypes but the distances of some of these points in the competitive class are longer than acceptable (D>ρ θi ). This causes rejecting a sufficient number of prototypes of the competitive class and the correct classification. In this case the classical kNN algorithm gives a wrong answer. Another possibility may occur when a correct and competitive class (classes) have an equally numerous representation among the prototypes of the shortest distance but the average distance of these points of the correct class is longer than that of the competitive one. In the situation where distances of a higher number of cases in the competitive class are longer than the acceptable kNN algorithm provides a wrong classification, whereas the proposed one correct. Sometimes there occurs a situation where the proposed for the adopted threshold rejects all the prototypes of both the correct class and the competitive one. In this situation we obtain a FRR error. Generally speaking, it is important to state that the kNN algorithm modified for the needs of biometrics is quite conservative and such a situation occurs relatively frequently. The classifier selects the answer 'I don't know' forcing the user to repeat the trial. The value of threshold ρ θi for particular people plays a very important role. A low value of the threshold increases the value of FRR error. Along with the increase in value of the threshold the modified kNN approaches in its work a classical kNN, which reduces FRR at the cost of higher FAR. The choice of the threshold value is a compromise between the desired flexibility of the system and its resistance. Thus, taking into consideration a threshold 0.7 we force a situation where it is slightly more often than every 20 person has to repeat an attempt at access to the system. However, this value secures only every 50th trial finishing with treating user X as user Y (FAR error 1.96%: 26 wrongly classified strides). It is important to note that in this case the number of wrong classifications is more than twice lower than in the case of the classical kNN algorithm. A certain anxiety is triggered by the results of the FAR error for intruders. At least 4% (for threshold 0.7 -9%) the acceptance of the unauthorized people is a result definitely unsatisfactory. In fact, the proposed classifier even in this case works better than a classical kNN (where the value of FAR error would be 100%); however, it indicates a need for seeking for further algorithm modifications, which will protect the biometric system from access of unauthorized people without a dramatic rise of FRR error of users.

Conclusion

The results obtained, basing on measurements of ground reaction forces generated during a walk conducted on an enormous, in comparison to other authors researches, data base independently from the applied classifier demonstrate a high potential of human gait as a biometric. The presented modification of the kNN algorithm for the needs of biometrics, to a considerable degree, fulfills its task. It allows for a correct classification of also those cases which have no chances with the classical version of kNN. Also a weak point of the proposed method was demonstrated: a relatively low resistance to attempts at access to resources on the part of unauthorized people. The whole encourages to continue works in this direction, which, for example more rigorous criteria of threshold ρ θi selection, should make a good method of classification still better.

Fig. 1 .

 1 Fig. 1. Components of GRF in: a) anterior/posterior, b) vertical, c) medial/lateral direction, during the support phase of the left lower limb

4

  Research Material and the Procedure of the Experiment 142 people (62 women and 80 men) took part in the research conducted in the Biaystok University of Technology. The participants in the research were 21.20 ± 1.14 years old, body weight 74.95 ± 17.0 kg and body height 174.26± 8.81cm.

Table 1 .

 1 Female and male characteristics

		No. of subjects	Age Body weight Body height
	Female	62 21.15±1.10 64.47±11.68 166.83±5.53
	Male	80 21.24±1.17 83.08±16.04 180.03±6.15
	5 Results and Disscusion
	Percentage of wrong classification for the classical kNN algorithm amounted
	FRR 2.19% (29 cycles) and FAR 4.00% (53 strides wrongly classified). So good an
	output results from a very big difference in the distances determined with DTW.

Table 2 .

 2 The results of the identification of the users and impostors using modified kNN

		modified kNN
	security level	users	impostors
	(thereshold)
		FRR FAR	FAR
	0.1	1.06% 3.25% 44.83%
	0.2	1.44% 3.17% 41.38%
	0.3	1.81% 2.87% 35.06%
	0.4	2.27% 2.79% 29.89%
	0.5	3.47% 2.57% 24.14%
	0.6	4.53% 2.34% 16.67%
	0.7	5.82% 1.96%	9.19%
	0.8	8.23% 1.66%	7.47%
	0.9	11.85% 1.06%	4.02%


Lionel Roques
email: lionel.roques@avignon.inra.fr

Marie-Anne Auger-Rozenberg

Alain Roques

Modelling the impact of an invasive insect, Megastigmus schimitscheki, via reaction-diffusion

Keywords: Insect dispersal, Reaction-diffusion model, Fickian diffusion, Fokker-Planck, Seed chalcid, Cedar

An exotic, specialist seed chalcid, Megastigmus schimitscheki, has been introduced along with its cedar host seeds from Turkey to southeastern France during the early 1990s. It is now expanding in plantations of Atlas Cedar (Cedrus atlantica). We propose a two-dimensional reactiondiffusion model to predict the expansion and impacts of this insect. Two main diffusion operators, of Fokker-Planck and Fickian types, are tested. We show that taking account of the dependence of the insects mobility with respect to spatial heterogeneity, and choosing the appropriate diffusion operator, are critical factors for obtaining good predictions.

Introduction

As a result of technological advances and of trade globalization, forest products are moved around the globe with an ever-increasing speed and frequency. This movement has precipitated a substantial increase in biological invasions by allowing organisms to pass the natural barriers that limit their dispersal (Liebhold et al., 1995). This is especially true for the insects infesting tree seeds, whose cryptic way of life does not usually allow a detection easier than using Xrays, the more as phytosanitary regulations are limited or non existent in many countries. In the West Palearctic, 8 of the 21 species recorded in the seed chalcid genus Megastigmus Dalman (Hymenoptera: Torymidae) are exotic invasive species introduced through seed trade [START_REF] Roques | Seed-infesting chalcids of the genus Megastigmus Dalman (Hymenoptera: Torymidae) native and introduced to Europe: taxonomy, host specificity and distribution[END_REF]. Modelling the expansion processes of such invasive species, and the related damages to forests and tree plantations, becomes a priority. [START_REF] Fabre | Competition between exotic and native insects for seed resources in trees of a mediterranean forest ecosystem[END_REF] just recorded a newly invasive seed chalcid, Megastigmus schimitscheki Novitzky (Hymenoptera: Torymidae) in southeastern France. The insect has probably been introduced from Asia Minor (Turkey) along with seeds of Cedar-of-Lebanon (Cedrus libani A. Rich.), during the period 1990-1994, in the region of Mont-Ventoux. The study of [START_REF] Fabre | Competition between exotic and native insects for seed resources in trees of a mediterranean forest ecosystem[END_REF], provides many useful data in order to build a model. In particular it is shown that the biological characteristics of M. schimitscheki in both native and invaded area are rather similar to these of other Megastigmus seed chalcids which have been extensively studied in Europe, North America and Japan (e.g. Kamijo, 1962;Hedlin et al., 1980;[START_REF] Roques | Les insectes ravageurs des cônes et graines de conifères en France[END_REF][START_REF] Roques | Seed-infesting chalcids of the genus Megastigmus Dalman (Hymenoptera: Torymidae) native and introduced to Europe: taxonomy, host specificity and distribution[END_REF]. Because of this recent introduction, M. schimitscheki did not recruit any predators in France yet. In addition, the insect is facing little competition for the exploitation of seed resources because the new host, Cedrus atlantica, is also introduced (but from North Africa) and its cones and seeds are still not colonized by the indigenous entomofauna. There is only one competitor, Megastigmus pinsapinis Hoffmeyer, another Megastigmus seed chalcid which has been introduced half a century ago from North Africa [START_REF] Roques | Seed-infesting chalcids of the genus Megastigmus Dalman (Hymenoptera: Torymidae) native and introduced to Europe: taxonomy, host specificity and distribution[END_REF]. However, M. pinsapinis seems dominated by M. schimitscheki [START_REF] Fabre | Competition between exotic and native insects for seed resources in trees of a mediterranean forest ecosystem[END_REF]. These patterns make the modelling of expansion and related damages easier than in the case of species introduced since a longer time and facing a large cortege of natural enemies and competitors. It allows in particular the use of single-species population dynamics models.

Reaction-diffusion models have proved themselves to be well adapted for describing some invasive species expansion in homogeneous and heterogeneous environments [START_REF] Skellam | Random dispersal in theoretical populations[END_REF][START_REF] Shigesada | Biological Invasions: Theory and Practice[END_REF], especially when one cannot follow the movement of each organism independently (for review, see [START_REF] Shigesada | Biological Invasions: Theory and Practice[END_REF][START_REF] Turchin | Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants[END_REF][START_REF] Murray | Mathematical biology[END_REF][START_REF] Okubo | Diffusion and Ecological Problems -Modern Perspectives[END_REF]. In the present case of M. schimitscheki, however, standard approaches may fail. Indeed, the species life-cycle is composed of two distinct stages: the ovo-larval stage (at least 2 years) and the adult stage (average duration of 10 days per individual). The dispersal occurs during the adult stage, which begins by progressive emergence of the adults from eggs of the previous years, and ends by the death of the adults. During this dispersal stage, laying -or, equivalently, attack-can occur, if the insects find free cedar seeds. To our knowledge, the methodology for coupling the ovo-larval stage time-discret model with a reaction-diffusion model for the adult stage, is not described in the existing literature. It requires to model adult progressive emergence and to relate the adult population density, obtained as the solution of the reaction-diffusion model, with cedar seeds attack rates.

The spatial heterogeneity considered in this paper is the presence or absence of cedar plantations. To deal with insect dispersion in such a patchy landscape, two main types of diffusion operators are found in the reaction-diffusion literature: the Fokker-Planck diffusion operator, and the Fickian diffusion operator. Fokker-Planck diffusion, which can be obtained as the macroscopic limit of uncorrelated random walks with spatially-dependent move lengths, is said to be more adapted than the Fickian one for modelling ecological processes in varying environments [START_REF] Turchin | Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants[END_REF]. However, Fickian diffusion is more adapted for solving theoretical problems, and both are widely used in the ecological modelling literature (see e.g. [START_REF] Shigesada | Biological Invasions: Theory and Practice[END_REF] for examples of Fickian diffusion, and Hannunen and Ekbom (2001) for an example of Fokker-Planck diffusion).

In Section 2, we propose a reaction-diffusion model for the adult stage, and a simple discrete model for the ovo-larval stage. These two models are coupled through a formula linking the solution of the reaction-diffusion model (adult population density) and the cedar seeds attack rates. In Section 3, Fokker-Planck and Fickian diffusion are discussed and further compared through numerical computations.

The model

Insect behaviour and simplifying hypotheses

Several biological traits facilitate the establishment of Megastigmus seed chalcids in exotic countries: (i) an ability to extend larval diapause in seeds for up to 5 years in most species, allowing them to bridge periods where seeds are unavailable because of masting phenomenon [START_REF] Turgeon | Insect fauna of coniferous seed cones: Diversity, Host plan interactions, and Management[END_REF]; (ii) parthenogenesis (the females can reproduce without males [START_REF] Roques | Les insectes ravageurs des cônes et graines de conifères en France[END_REF][START_REF] Roques | Seed-infesting chalcids of the genus Megastigmus Dalman (Hymenoptera: Torymidae) native and introduced to Europe: taxonomy, host specificity and distribution[END_REF]. A simplified life cycle of M. schimitscheki is given in Fig. 1. The data obtained by [START_REF] Fabre | Competition between exotic and native insects for seed resources in trees of a mediterranean forest ecosystem[END_REF] and [START_REF] Boivin | Potentialités invasives et génétiques des populations d'insectes ravageurs de graines de Conifères : comparaison de deux espèces du genre Metastigmus[END_REF] revealed about the same proportion of females and males at emergence (mean male/female ratio: 0.99, standard deviation: 0.28). This allows us to neglect the fecundation processes, and to focus on female dispersion only.

Female emergence occurs over a 6-11-day period, with a quite constant rate (Fabre et al., 2003, Fig. 2). Similarly as in the related species, Megastigmus spermotrophus, the mean individual life expectancy is assumed to last 10 days (Jarry et al., 1997).

The chalcid M. schimitscheki is considered as specifically related to seeds of cedars (Cedrus spp.). Furthermore, behavioural observations using videorecording showed that adults having landed on isolated host trees bearing non or only a few cones usually fly up, in a slow vertical movement, to the tree top, and then are blown away from the tree (A Roques, unpublished observations). Visual and olfactary attractiveness of cedar cones have also been observed [START_REF] Roques | Response of Megastigmus spermotrophus to colour traps baited with volatile constituents of Pseudotsuga menziessi seed cones[END_REF]Luik et al., 1999). Thus, we make the assumption that, inside the cedar plantations, the adults remain at low altitude, while outside the plantations, they fly at a higher altitude, and are therefore subject to stronger winds and air turbulences [START_REF] Nathan | Mechanistic models for tree seed dispersal by wind in dense forests and open landscapes[END_REF].

Neglecting environment topography, and the forest and plantations other than cedar, we consider that the environment is made of two components: cedar plantations and other regions.

2.2

Model for the adult stage Note: In this section, insect emergence is assumed to occur during year n.

Reaction-diffusion models are known to work properly in some situations of redistribution of small organisms. In the model proposed below, the reaction part of the equation is not totally standard. Indeed, it does not have to take the reproduction into account, and thus differs from the classical logistic nonlinearity u(ε -νu). However, it has to take account of the insects progressive emergence and of the individuals mortality. As we mentioned in the introduction of this paper, two types of diffusion operators are considered, leading to two models. The model with a Fokker-Planck diffusion:

∂u ∂t (t, x) = ∇ 2 (D(x)u) - u ν + f (t, x), (2.1)
and the model with Fickian diffusion:

∂u ∂t (t, x) = ∇(D(x)∇u) - u ν + f (t, x). (2.2)
In order to assess the usefulness of considering such models with a spacedependent coefficient D(x), we also introduce the model with constant diffusion coefficient

D 0 : ∂u ∂t (t, x) = D∇ 2 u - u ν + f (t, x). (2.3) In each model, t ∈ [0, N d ].
The initial time t = 0 corresponds to the beginning of emergence, and N d is an arbitrary time, sufficiently large such that living adults cannot be found at t = N d . The unknown u(t, x) is the female adults population density at time t and spatial coordinate x = (x, y). These equations are considered in a bounded domain Ω ⊂ R 2 , with lethal boundary conditions u = 0 on ∂Ω. The domain Ω is chosen large enough so that taking lethal boundaries is realistic, since in fact, almost no individual reaches the boundary. The coefficient [START_REF] Turchin | Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants[END_REF]. The idea of considering such spatiallydependent motility and diffusivity coefficients was motivated by the biological observations ( § 2.1) regarding the behaviour of M. schimitscheki. Indeed, if the insects increase their flight altitude as they leave a cedar plantation, they become subject to stronger winds and turbulences [START_REF] Nathan | Mechanistic models for tree seed dispersal by wind in dense forests and open landscapes[END_REF], and their mobility is thus increased. Though the wind effects are usually modelled by adding an advective term of the form v(t, x)•∇u to the equation, it appears here to be more judicious to incorporate these effects into the "random" dispersion term, since the local wind cannot be measured in practice, and may take very different directions, in function of the local topography. See Section 3 for a comparison between models (2.1), (2.2) and (2.3).

D(x) in model (2.1) is called motility. In model (2.2), D(x) is called diffusivity)
The mortality rate is assumed to be constant in time. This is modelled by the term -u ν , which means that, at each time unit, a proportion 1/ν of the individuals dies. A straightforward computation leads to ν = life expectancy.

Note:

The dispersion operators considered in these models are conservative. However, because of the mortality term, -u/ν, and of the lethal boundary conditions, the population density converges to 0 as t increases, and is almost 0 at t = N d , which is consistent with the definition of N d .

The progressive emergence of the individuals occurs during a time period [0, N em ]. It is modelled by the term f (t, x). Defining the scalar field m(x) as the density, before emergence, of female adults that will emerge during year n, at the position x ∈ R 2 , we have

f (x, t) = m(x) N em for t ∈ [0, N em ], and f (x, t) = 0 for t > N em .
(2.4)

The underlying assumption in (2.4) is that the emergence rate is constant (see § 2.1). The function m(x) is computed as follows. Let n p be the number of cedar plantations in Ω. We denote these plantations by P i , for i = 1 . . . n p , and by A(P i ) their respective areas. The function m(x) is given by the formula:

m(x) = U n (P i ) A(P i ) .
(2.5)

Here U n (P i ) corresponds to the number of female adults that will emerge during year n in plantation P i , and is given by the ovo-larval stage model (see § 2.4).

Computation of the attack rates

To our knowledge, no standard method exists for computing attack rates (or other type of impacts), from the population density solution of a reactiondiffusion model. We derived a formula linking the attack rates in each plantation with the solution u of the model and some biological parameters. Let P i be a given plantation, and E i be the expectancy of the number of new eggs laid in this plantation. Let U n = n p j=1 U n (P j ), be the total number of female adults which emerge during year n. For large values of U n (for M. schimitscheki, we generally deal with values of U n above 10 4 ), and assuming that the individuals are independent, E i can be computed in the following way:

E i = εT i U n , (2.6)
where ε is the mean laying frequency (number of eggs per unit of time) and T i is the mean time spent in plantation P i . The number of individuals in plantation P i , at time t, is equal to Pi u(t, x)dx. Thus, the mean time spent by each individual in P i is:

T i = N d 0 P i u(t, x)dxdt U n .
(2.7)

Combining (2.6) and (2.7), we get:

E i = ε N d 0 Pi u(t, x)dxdt.
(2.8)

Thus, the estimated attack rate A n (P i ), of year n, in the plantation P i , at the end of the adult stage, is:

A n (P i ) = ε F n (P i ) N d 0 P i u(t, x)dtdx, (2.9)
where F n (P i ) corresponds to the fructification of year n (number of seeds) in plantation P i . Note that, in practice, F n (P i ) is generally measured in the course of year n + 1. Indeed, the attack occurs one year before the harvest takes place.

Model for the ovo-larval stage

Modelling this part of the insect life cycle does not require sophisticated mathematical tools. However, we briefly present a model here, for the sake of completeness. For each year k, the fructifications (number of seeds) and attack rates (percentage of seeds infested by M. schimitscheki ) in each plantation P i are denoted by F k (P i ) and A k (P i ), respectively. The number U k (P i ) of female adults that will emerge, during year n, in the plantation P i is defined recursively by:

U n (P i ) = 5 k=2 σ n-k A n-k (P i )F n-k (P i ) 2 , ( 2.10)
where σ n-k is a factor corresponding to the ratio of individuals laid during year n -k, which emerge during year n. This factor takes account of the prolonged diapause rates (see § 2.1), and the mortality during the ovo-larval stage. In accordance with empirical observations, we assumed here that the ratio female/male was 1, hence the divisor 2.

3 Comparison of models (2.1) and (2.2), numerical computations

As recalled in the introduction of this paper, the Fokker-Planck diffusion term ∇ 2 (D(x)u) of model (2.1), seems to be better-adapted to the modelling of insects dispersal than the Fickian diffusion term ∇(D(x)∇u) of model (2.2) which is, on the other hand, well adapted for physical problems such as the description of the heat propagation in an inhomogeneous medium. Assume that we are in a binary environment, and that the coefficient D(x) is 2 times differentiable and takes the constant values a > 0 inside the plantations, and b > 0 outside the plantations, with b > a from the observations of § 2.1. We have

∇ 2 (D(x)u) = ∇(D(x)∇u + u∇D(x)), (3.11)
thus, in the model with Fokker-Planck diffusion, an additional convection term ∇(u∇D(x)) appears compared to the model with Fickian diffusion. In regions where D(x) is constant, this term vanishes, and the dispersion operators are therefore equivalent. However, at the interface between plantations and outside the plantations, with this additional term, the individuals in the first model with Fokker-Planck diffusion, are more likely to go from the regions where the motility D(x) is high to the regions where D(x) is low. Thus, the first model should predict higher attack rates than the second one, at least in the regions of emergence. Note that, when D is constant in Ω, models (2.1) and (2.2) are both equivalent to model (2.3).

In order to get a better idea of the differences between the models (2.1), (2.2) and (2.3), we carried out numerical simulations of the adult expansion and related attack rates, on an hypothetical landscape distribution. Most of the parameters values were deduced from previous empirical studies.

The landscape distribution that we considered is presented in Fig. 2. We considered a square domain Ω, of sides 100 km, and including four cedar plantations of same circular shape, with a radius R = 0.5 km. This lead to areas of 78 ha, which are realistic in the context of southeastern France cedar plantations. We labelled these plantation P 0 , P 1 , P 2 and P 3 , and we assumed that the number of emerging females U n (P i ), i = 0 . . . 3 was known. Namely, we took U n (P 0 ) = 10 5 , and U n (P i ) = 0 for i = 1 . . . 3. Thus P 0 was assumed to be the "source" point of the insect expansion. Plantation P 1 and P 2 were placed at 5 km and 15 km, respectively, to the East of P 0 . Plantation P 3 was placed at 15km to the South of P 0 . The fructification F n (P i ), i = 0 . . . 3, was assumed to be the same in each plantation, and equal to 10 6 seeds, which is reasonable for plantations of this size (data obtained from the French National Forestry Office). The the initial ratio (emerging males+females)/(available seeds) is therefore of 0.2 in P 0 .

The mortality rate, 1/ν, was directly computed, using a life expectancy, ν of 10 days (Jarry et al., 1997). The mean laying rate ε = 1 was obtained from previously conduced experiments on M. schimitscheki [START_REF] Boivin | Potentialités invasives et génétiques des populations d'insectes ravageurs de graines de Conifères : comparaison de deux espèces du genre Metastigmus[END_REF]. The length of the emergence period was set to N em = 9 days (cf. § 2.1), and N d was set to 40 days. In models (2.1) and (2.2), the only parameter which was not known precisely was the function D(x) (or the constant D 0 for model (2.3)). We assumed that D(x) was constant equal to a value a > 0 inside the plantations, and was also constant but equal to a certain value b > a, at a distance 2R = 1 km of the plantations centres. In the rings comprised between a distance R and 2R from the plantations centres, we assumed that D(x) was radially increasing, from the centre of the closest plantation. More precisely, setting r(x) = the distance between x and the centre of the closest plantation, we have:

   D(x) = a if r(x) ≤ R, D(x) = a + (b-a)(r(x)-R) 3 R 5 (6r(x) 2 -27r(x)R + 31R 2 ) if R ≤ r(x) ≤ 2R, D(x) = b if r(x) ≥ 2R.
(3.12) This function is depicted in Fig. 3.

In models (2.1) and (2.2), we assumed that the value of a was 0.1 km 2 /day, corresponding to a low rate of movement inside the plantations (see [START_REF] Shigesada | Biological Invasions: Theory and Practice[END_REF] for some observed values of D, for different species). The value of b was set either to 1 km 2 /day or to 5 km 2 /day, corresponding respectively to intermediate, and high rates of movement. In the case of model ( 2.3), we set either D 0 = 0.1 km 2 /day, D 0 = 1 km 2 /day or D 0 = 5 km 2 /day.

We solved the reaction-diffusion models using a finite elements method, using the Comsol c software.

Results

For each model, and each parameter value, we numerically obtained a population density u(t, x), for t ∈ [0, N d ] and x ∈ Ω. The corresponding attack rates A n (P i ), in plantations P i , i = 0 . . . 3 were computed thanks to formula (2.9).

In this section, we focus our analysis on the attack rates rather than on the population density, since, in practice, only the attack rates are measured. Fig. 4 gives an overview of the population densities solutions of our models, in P 0 and P 1 , for particular values of the parameters. For model (2.1), when a = 0.1 and b = 1, the attack rates in plantations P 0 , P 1 , P 2 and P 3 at the end of the invasion are 19.78%, 1.18%, 0.02% and 0.02% respectively. When b = 5, these attack rates are 26.48%, 3.80%, 0.44% and 0.51% respectively (Fig. 5,(a)). The attack rates are larger for b = 5, even in the source plantation P 0 . This phenomenon is probably the consequence of an higher convection coefficient, oriented inside the plantations, for large values of b. In both cases, the observed attack rates are consistent with what could be expected form a biological viewpoint. Note that the attack rates in plantations P 2 are only slightly smaller than in P 3 , which advocates for a low barrier effect played by the plantation P 1 .

For model (2.2), when a = 0.1 and b = 1, the attack rates are 7.67%, 0.21%, 0.00% and 0.00%, respectively in plantations P i , i = 0 . . . 3. For b = 5, we obtained 5.31%, 0.15%, 0.02% and 0.02% (Fig. 5,(b)). Thus, in the case of model (2.2), when b is increased, the plantations P 0 and P 1 have lower attack rates, and the plantations P 2 and P 3 , which are farther from the invasion source point, have slightly higher attack rates. In both cases, however, the attack rates are very small compared to model (2.1), and do not seem realistic from a biological viewpoint, basing ourselves on preliminary experiments conduced in the region of Mont-Ventoux. Indeed, these attack rates are too low to sustain an invasion (compare with the initial ratio (emerging males+females)/(available seeds) of 0.2 in P 0 ). Lastly, for model (2.3), when D 0 = 0.1, the attack rates are 15.21%, 0.04%, 0.00% and 0.00% respectively in P i , i = 0 . . . 3. When D 0 = 1, it gives 2.87%, 0.24%, 0.01% and 0.01%, and for D 0 = 5, 0.77%, 0.17%, 0.02% and 0.02% (Fig. 5,(c)). In all cases, the attack rates given by this model are smaller than those given by model (2.1). When D 0 = 0.1, the attack rate in P 0 remains at a realistic level, but the other attack rates are very low, compared to what could be expected. For higher D values, the attack rate in the source point P 0 either, is too low to sustain an invasion.

The model (2.1), seems better-adapted to our problem than the two other models. We also tested on this model the influence of the length of the emergence period on the attack rates. With b = 5, we obtained, taking N em = 1, attacks rates of 25.73%, 3.72%, 0.44% and 0.51%, respectively in plantations P i , i = 0 . . . 3. Thus, in this case corresponding to a massive emergence of all the individuals in one day, the attack rates are very slightly decreased compared to the case of an emergence over a nine-day period.

Discussion

In this work, we have proposed models to describe the expansion of an invasive insect, and its impact on cedar plantations. The life-cycle of the insect consisted in two distinct stages: ovo-larval and adult stages. Modelling the ovo-larval stage of the insect did not require sophisticated mathematical tools. On the other hand, the insect expansion, occurring during the adult stage was modelled with a two-dimensional reaction-diffusion equation. Though reaction-diffusion models are very common in the literature related to insect dispersal (see the books [START_REF] Shigesada | Biological Invasions: Theory and Practice[END_REF][START_REF] Turchin | Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants[END_REF][START_REF] Murray | Mathematical biology[END_REF][START_REF] Okubo | Diffusion and Ecological Problems -Modern Perspectives[END_REF], for review), we raised some questions which, to our knowledge, are not classical.

The first point was to relate the models for the ovo-larval stage and adult stage. To do so, we proposed a simple formula, in § 2.3, linking the adult population density, solution of a reaction-diffusion model, with a seed attack rate, corresponding to the new laid eggs at the end of the seasonal invasion (Fig. 1). The idea behind this formula is that the attack rate in each plantation is proportional to the mean time spent by the individuals inside this plantation.

The second point was to chose the appropriate diffusion operator. We considered the three most common diffusion operators found in the reaction-diffusion literature. The corresponding models were (2.1), with Fokker-Planck diffusion, (2.2) with Fickian diffusion, and the model with constant diffusion coefficient D. The wind effects were not modelled through an additional advection term, since, at the scale considered here, winds may vary greatly in time and space, and cannot therefore be measured precisely. We chose rather to incorporate them in the random diffusion term. We assumed, in models (2.1) and (2.2), that the mobility was smaller inside the plantations. This assumption was consistent with previous empirical observations, since Megastigmus fly at lower altitudes inside the plantations, and are therefore subject to lower winds.

Writing the Fokker-Planck diffusion operator of model (2.1) in a different form, we observed that the main difference with model (2.2), and also with model (2.3), was a convection term, pushing the individuals situated at the interface of the plantations into these plantations. Numerical computations confirmed this observation. With the same coefficient D(x), model (2.1) predicted higher attack rates than model (2.2). The model (2.3), with constant coefficient D 0 , predicted very low attack rates, except near the invasion source point for small D 0 values. The comparison with model (2.1) showed that taking account of the effects of spatial heterogeneities on the individuals mobility, could have an important impact on the predicted attack rates. From our point of view, and basing ourselves on some preliminary observations made at Mont-Ventoux, model (2.1) gave more realistic results than the two others, for the species considered here. The model (2.2) tended to "homogenize" the population density, has it does for the solution of the heat equation. Indeed, as soon as the emergence period ends, almost no individual remains in the source plantation (see Fig. 4,(b)). To correct this, a drift term, oriented towards the region where D is smaller, should be added. This would lead to models close to (2.1). We would like to recall, however, that the "self-adjoint" nature of the Fickian diffusion operator considered in model (2.2) allows to obtain important theoretical results (see [START_REF] Cantrell | Spatial Ecology via Reaction-Diffusion Equations[END_REF]Kinezaki et al., 2003;Berestycki et al., 2005a,b) and can be more realistic in other situations.

In order to fit the unknown coefficients a and b, and to assess the validity of the model (2.1), we are planing to obtain data from the region of Mont-Ventoux. The location of the cedar plantations in Southeastern France will be obtained from the French National Forestry Office (FNFO). On each plantation, the FNFO will give estimations of the annual seed production in each plantation. Then, standardized seeds samplings will be carried out. Immediately after being collected the seeds will be X-rayed in order to calculate the proportion of filled, empty and chalcid-infested seeds. However, M. schimitscheki cannot be separated from M. pinsapinis through observations of larvae on Xray plates. Therefore, the infested seeds will be placed in emergence boxes until adult emergence in order to know the proportion of each species, and to estimate the specific damage that can be attributed to M. schimitscheki.

In this paper, we have explicitly modelled adults progressive emergence. We have seen in § 3.1 that, when the length of the emergence was shortened, the model (2.1) did not give dramatically different results. It is not surprising since no interaction between the individuals have been assumed. However, modelling this progressive emergence remains a very important question. Indeed, the domination of the competing species M. pinsapinis by M. schimitscheki seems to be linked to the earlier emergence of M. schimitscheki [START_REF] Fabre | Competition between exotic and native insects for seed resources in trees of a mediterranean forest ecosystem[END_REF]. In a future work, we plan to build and analyse the competition between these two species through a system of reaction-diffusion. Much of the tools used in this paper could surely be adapted to this case. Taking account of the shape of the emergence function will probably be of crucial interest.

Lastly, we would like to emphasize that other types of models exist and are efficient in practice for modelling such dispersion phenomena in heteroge-neous environments. For instance, individual based models (see e.g. Kareiva and Shigesada, 1983;[START_REF] Marsh | The form and consequences of random walk movement models[END_REF][START_REF] Gross | Olfactory responses of seed wasps Megastigmus pinus Parfitt and Megastigmus rafni Hoffmeyer (Hymenoptera Torymidae) to host tree odours and some monoterpenes[END_REF], where each individual interacts with its environment, and simulations are run until a general behaviour of the population can be observed. Integro-difference or dispersal kernel models can also be useful especially when long distance dispersal is a key factor of the expansion [START_REF] Mollison | Spatial contact models for ecological and epidemic spread[END_REF]Klein et al., 2006). However, our aim in this paper was precisely to compare some dispersion operators used in practice in reaction-diffusion models, and to give a clear approach of the modelling of insects expansion and impacts in the framework of reaction-diffusion equations.

The research was partly funded by a grant from the French Ministry of Agriculture "DERF-INRA N61.45.19/01-Caractérisation potentialités invasives d'un nouveau ravageur des graines de cèdre récemment introduit au Mont Ventoux". The research has also been partly funded by the European Commission within the FP 6 Integrated Project "ALARM" (GOCE-CT-2003-506675). 2.2). The boundary of the domain Ω, whose center is assumed to be located in P 0 , is not represented here.

Figure captions

Fig. 3: The function D(x), defined by (3.12), near a circular plantation of radius R = 0.5, and centered in the point (0, 0). In this plot, it is assumed that a = 0.1 and b = 1.

Fig. 1 :

 1 Fig. 1: Simplified life-cycle of M. schimitscheki.

Fig. 2 :

 2 Fig.2: Hypothetical plantations distribution used for the numerical computations. The plantations are represented by plain circles. The dashed circles delimit the regions where D(x) is not constant in models (2.1) and (2.2). The boundary of the domain Ω, whose center is assumed to be located in P 0 , is not represented here.

Fig. 4 :

 4 Fig. 4: Mean population densities in plantations P 0 (solid lines) and P 1 (dashed lines), during the time-span [0, N d ]. (a): Model (2.1), with a = 0.1 and b = 5. (b): Model (2.2), with a = 0.1 and b = 5. (c): Model (2.3), with D 0 = 1.

Fig. 5 :

 5 Fig. 5: Bar charts presenting the predicted attack rates, in each plantation. (a): Model (2.1); the black bars correspond to the case b = 1, and the white bars correspond to b = 5. (b): Model (2.2); the black bars correspond to the case b = 1, and the white bars correspond to b = 5. (c): Model (2.3); the black bars correspond to the case D 0 = 0.1, the gray bars correspond to D 0 = 1, and the white bars to D 0 = 5.

  Figure 2:

Figure 4 :

 4 Figure 4:

Figure 5 :

 5 Figure 5:

Acknowledgements

Ovo-larval stage, inside the cedar seeds.

Usual duration: 2 years. For a part of larvae: 3 to 5 years ("prolonged diapause").

Seed of C. atlantica infested by a Megastigmus larva (X-ray pic.)

Year n, May-Year n, Jun M. schimitscheki adult.

Some insects find trees with non-infested seeds. Some insects do not find non-infested seeds, and die.

Year n, May-Year n, Jun

Year n, May-Year n, Jun

Each female lays approximately 1 egg/day.

Year n, May-Year n, Jun

The adults die. The eggs become larvae.

Emergence period: 6-11 days. The emergence rate is almost constant. Adults mean life expectancy: 10 days.


Dariusz Frejlichowski
email: dfrejlichowski@wi.zut.edu.pl

Piotr Czapiewski
email: pczapiewski@wi.zut.edu.pl

An Application of the Curvature Scale Space Shape Descriptor for Forensic Human Identification Based on Orthopantomograms

Keywords: dental images, forensics, shape descriptors, Curvature Scale Space

This paper is devoted to the problem of forensic human identification based on the application of dental biometrics. For this purpose digital orthopantomogram images are used. The complete scheme of the approach is described. It starts with the image enhancement by means of the Laplacian pyramid decomposition. Secondly, the segmentation of an image into sub-images covering single teeth is performed. The process is based on determining the line between the upper and lower jaws by means of integral projections and, later, image intensity and types of teeth. For each obtained sub-image, covering a single tooth, its shape is extracted and the outer contour of this shape is traced. Thus, a shape descriptor can be applied in order to represent an object for the later identification. In this paper, the Curvature Scale Space algorithm is applied and investigated as one of the elements of the biometric identification approach based on digital orthopantomogram images.

Introduction

Dental records are widely used in forensic human identification as one of the most efficient biometric modalities in this application [START_REF] Bowers | Forensic Dental Evidence[END_REF]. They are often used in the process of establishing the identity of a person, together with DNA and fingerprints. The popularity of dental biometrics in judicial proceedings comes from their property of discrimination and the robustness to decomposition [START_REF] Lee | The Diversity of Dental Patterns in Orthopantomography and its Significance in Human Identification[END_REF], where for example automatic face recognition would fail [START_REF] Spaun | Face Recognition in Forensic Science[END_REF][START_REF] Shchegoleva | Facial Surface Reconstruction in 3D Format[END_REF]. It has led to several attempts for the development of an identification system, such as WinID that compares dental records previously codified by an expert. The second example, which is closer to the approach assumed in the works described in this paper, is the Automated Dental Identification System (ADIS). One of the stages in this system is the automatic extraction of dental biometrics [START_REF] Fahmy | Towards an Automated Dental Identification System (ADIS)[END_REF]. Whereas WinID applies dental works as a basis for identification, ADIS utilizes teeth shapes extracted from radiograms for this purpose. The ADIS system is an important tool used in the work of the Dental Task Force (DTF), a special service created by the Federal Bureau of Investigation (FBI) [START_REF] Nassar | A Prototype Automated Dental Identification System (ADIS)[END_REF].

The analysis of radiographic images is extremely significant for various medical and forensic purposes, and hence many general methods related to radiograms have been described [START_REF] Gut | Edge-Based Robust Image Registration for Incomplete and Partly Erroneous Data[END_REF][START_REF] Bator | Finding Regions of Interest for Cancerous Masses Enhanced by Elimination of Linear Structures and Considerations on Detection Correctness Measures in Mammography[END_REF][START_REF] Kukolowicz | Evaluation of Set-up Deviations During the Irradiation of Patients Suffering from Breast Cancer Treated With Two Different Techniques[END_REF]. However, the specifics of the dental image, and especially of orthopantomograms, must be taken into consideration in order to properly discern between teeth and represent their shape.

The automatic system for the identification of persons based on digital orthopantomograms is composed of image enhancement and segmentation, and feature extraction and recognition. In the first step one can apply, for example, the wavelet-transform [START_REF] Lu | Contrast Enhancement of Medical Images Using Multiscale Edge Representation[END_REF][START_REF] Dippel | Multiscale Contrast Ehnahncement for Radiographies: Laplacian Pyramid Versus Fast Wavelet Transform[END_REF] or Laplacian pyramid [START_REF] Zhou | A Content-Based System for Human Identification Based on Bitewing Dental X-ray Images[END_REF][START_REF] Frejlichowski | Application of the Laplacian Pyramid Decomposition to the Enhancement of Digital Dental Radiographic Images for the Automatic Person Identification[END_REF]. During the segmentation process the integral projections [START_REF] Jain | Matching of Dental X-ray Images for Human Identification[END_REF][START_REF] Frejlichowski | Automatic Segmentation of Digital Orthopantomograms for Forensic Human Identification[END_REF] or active contour models [START_REF] Zhou | A Content-Based System for Human Identification Based on Bitewing Dental X-ray Images[END_REF] have been applied so far. The extraction of particular teeth shapes is performed using active shapes [START_REF] Chen | Automatic Forensic Dental Identification[END_REF], line scanning [START_REF] Jain | Matching of Dental X-ray Images for Human Identification[END_REF] or watersheds [START_REF] Frejlichowski | Extraction of Teeth Shapes from Orthopantomograms for Forensic Human Identification[END_REF]. The first two algorithms were developed specifically for the intraoral images. Hence, they fail when applied in pantomograms. This problem results from frequent occlusions. The third mentioned method performs better in this case. In the last stage any of the shape description and recognition approaches can be applied.

In this paper the whole algorithm is described. Firstly, the methods for image enhancement are presented. Next, the segmentation and localisation of the object of interest -single teeth -is provided. Both methods are discussed in the second section. Later, in the third section, the algorithm for the representation of the extracted contours is described. The fourth section discusses the experimental conditions and results. Finally, the last section concludes the paper.

Image Enhancement, Segmentation and the Extraction of Teeth Shapes

The image enhancement applied for the pantomograms starts with the Laplacian decomposition, proposed in [START_REF] Burt | The Laplacian Pyramid as a Compact Image Code[END_REF]. The Laplacian pyramid's layers are derived through the subtraction of consecutive layers of a Gaussian pyramid. This process is represented by means of the following formula:

X k =↓ ( Xk-1 ), L k = X k-1 -↑ (X k ), (1)
where X 0 is the original image, ↓ (X) and ↑ (X) represent the process of downsampling and upsampling the image by a factor of 2, Xk is a low-pass filtered image of X k , and L k is the successive layer of the Laplacian pyramid. The lowpass filtering is performed using the Gaussian filter.

For the Laplacian pyramid, a contrast equalisation function is used (as described in [START_REF] Vuylsteke | Image Processing in Computed Radiography[END_REF]):

f (x) = a( x |x| )|x| p , ( 2
)
where "x is normalized to the range [-1, 1] and the factor a is needed for rescaling the resulting image to the original dynamic range" [START_REF] Vuylsteke | Image Processing in Computed Radiography[END_REF]. The exponent p controls the degree of non-linearity. According to further modifications provided by [START_REF] Stahl | Noise-Resistant Weak-Structure Enhancement For Digital Radiography[END_REF], a new representation of the Laplacian pyramid, applied in our approach, can be formulated as:

{ r(x) = G • x • (1 -|x| M ) p + x, if |x| ≤ M r(x) = x, elsewhere , ( 3
)
where M is the upper limit for linear enhancement and G is the constant gain.

In our previous experiments [START_REF] Frejlichowski | Application of the Laplacian Pyramid Decomposition to the Enhancement of Digital Dental Radiographic Images for the Automatic Person Identification[END_REF][START_REF] Frejlichowski | Automatic Segmentation of Digital Orthopantomograms for Forensic Human Identification[END_REF][START_REF] Frejlichowski | Extraction of Teeth Shapes from Orthopantomograms for Forensic Human Identification[END_REF] the best results were obtained when using the following combination of methods: 1) the averaging of the two Laplacian pyramid layers next to last, 2) unsharp filter on the second layer, 3) contrast enhancement, based on the eq. 2-3. This approach was used at the stage of image pre-processing.

The second stage consists of the segmentation of an image. Firstly, the upper and lower jaws are separated by means of integral projections. Usually, the initial starting point for this process is established by the user. However, our goal is to fully automate this process. Hence, in the method described in this paper, this point is selected as the horizontal integral projection around the centre of the image with the lowest value. Taking into account the characteristics of a pantomogram (the separating line is not horizontal), only some of the pixels (20% of the image width), which lay closest to the frontal teeth gap selected earlier, are used for the derivation of projections. The area between the necks of two adjacent teeth is easy to locate on the pantomogram. Thus, the aforementioned line is applied for the location of the position of the neck for each tooth. Later, points lying on splines representing gaps between the necks of two adjacent teeth are selected. For the upper and lower jaw, an array of values is obtained, comprised of the intensity values of the points lying on the splines going through their respective dental pulps. The sharp spikes visible on the plot representing these values correspond to the dark regions surrounded by light intensities, what indicates a gap between teeth necks. For the matching process only small subsets (with the pre-assumed number of elements) of the values in the array are selected.

For the indication of a spike value corresponding to the gap between the necks of teeth for a processed tooth, lying on the position p c , the Bayesian probability P (x i , p c ) calculated for each point on the curve (x i ) is used:

P (x i , p c ) = I(x i )G(x i , p c )D(x i , p c ), (4)
where: I -the intensity of the range filtered and inverted original image in the point x i , G -a discrete Gaussian function with expected value equal to the horizontal position of the last detected gap, displaced left or right (depending on the current search direction) by the amount of pixels that is equal to the average width of the tooth at position p c , D -a function introduced to limit the number of pixels analysed in each iteration. It is equal to 1 for points, for which the horizontal distance from the previous detected gap is between 75% and 175% of the expected width of a tooth on position p c , with regard to the current search direction, and equal to 0 elsewhere. At each iteration, the current maximal value obtained using the above formula is added to the list of gap positions and it is used as a starting point for the next iteration. The iterations end at the border of an image or when the number of gaps for either the upper or lower jaw is equal to 8. Usually, the obtained gap locations are reliable indicators for the space between the teeth. However, sometimes the vertical line is insufficient for the distinction between two occluding teeth. Thus, additionally another point between them is localised by means of a greedy algorithm, that iteratively moves one pixel towards the radiogram's top or bottom, selects the pixel in horizontal neighbourhood with the highest intensity on the inverted, range-filtered image and uses it as the basis for the next iteration. The number of iterations is equal to the half length of an average tooth on a pantomogram. The last obtained point is taken as the second point applied for the separation. Using it, as well as the previously derived point, the segmentation line is established.

Later, the areas lying below the teeth roots are removed by means of the curve separating jaws. This line is shifted vertically until such an alignment is obtained, where the sum of pixels through which it passes is smaller than the surrounding result. That indicates the area between the teeth line and the cheekbone.

An exemplary pantomogram enhanced and segmented by means of the described approach is provided in Fig. 1. After the segmentation of the image, we can extract the shapes of particular teeth. Firstly, in order to lower the influence of noise and enlarge the areas with similar intensity, the opening operation on the image is performed. Later, the entropy filtering is applied. It detects edges of areas with similar colours. Then, the watershed algorithm is used in order to divide the image into small segments. The dimensions of the obtained areas depends on the size of a structuring element used in the morphological opening -the larger the element, the larger size of the segments.

For each segment the following features are obtained: the centroid, the normalised mean intensity, and the normalised vertical distance from the centroid to the curve separating the jaws. Next, the distance between obtained centroids is derived, and 50 of the closest segments to a processed segment are selected for the calculation of the average intensity. The distinction value of the i-th segment is derived using the equation:

D(i) = N ∑ j=1 max( Ī(i) -Ī(j), 0), (5)
where: Ī(i) and Ī(j) are the average intensities of i-th and j-th segments. These values help in determining whether a segment is brighter than its surrounding segments or not. Later, a fitness function is calculated in order to select the segments belonging to the tooth. The values of this function depend on the type of the tooth. For the incissors only the distinction function and vertical distance from the jaws separating curve are used, with the weights of 0.7 and 0.3, respectively. For other types of teeth the average intensity is also included, with the weights of 0.4 (distinction function), 0.4 (average intensity) and 0.2 (vertical distance from the jaws separating curve).

When all fitness functions for segments are obtained, the process of selection of those belonging to a tooth begins. Basing on the performed experiments we could assume that a fitness function value has to be larger than 0.4 for incissors, and 0.5 for the canine and first premolar. For other types of teeth the threshold is established as the average intensity multiplied by 0.8. The selected regions are morphologically dilated, what results in the removal of borders that separate them. Finally, the contour of a tooth can be traced and extracted (see Fig. 2).

Representation of the Teeth Contours Using the Curvature Scale Space Transform

The extracted contours of teeth are described using the Curvature Scale Space (CSS) transform [START_REF] Mokhtarian | Silhouette-Based Occluded Object Recognition Through Curvature Scale Space[END_REF][START_REF] Roh | Robust Player Gesture Spotting and Recognition in Low-Resolution Sports Video[END_REF]. The CSS representation is calculated using the convolution of a path-based parametric representation of a curve extracted from the boundary of processed shape with a Gaussian function of increasing variance σ 2 . The zeros of curvatures of the convolved curves are extracted and combined in a scale space representation. Those values are derived during the evolution of the planar curve changed by the expanding Gaussian function.

The closed planar curve r is represented parametrically for axes x and y using the normalized arc length parameter u:

r(u) = {x(u), y(u)|u ∈ [0, 1]}. (6)
The evolved curve is represented as Γ σ :

Γ σ (u) = {χ(u, σ), ψ(u, σ)}, (7)
where:

χ(u, σ) = x(u) ⊗ g(u, σ), ψ(u, σ) = y(u) ⊗ g(u, σ), (8)
and: ⊗ -the convolution operator, g -Gaussian function of width σ, calculated using the formula:

g(u, σ) = 1 σ √ 2π e -u 2 /2σ 2 . ( 9
)
The curvature κ of Γ can be formulated as:

κ(u, σ) = χ u (u, σ) -ψ uu (u, σ) -χ uu (u, σ) -ψ u (u, σ) (χ u (u, σ) 2 + ψ u (u, σ) 2 ) 3/2 , ( 10
)
where:

χ u (u, σ) = ∂ ∂u (x(u) ⊗ g(u, σ)) = x(u) ⊗ g u (u, σ), (11)
χ uu (u, σ) = ∂ 2 ∂ 2 u (x(u) ⊗ g(u, σ)) = x(u) ⊗ g uu (u, σ), ( 12
)
ψ u (u, σ) = y(u) ⊗ g u (u, σ), ( 13
)
ψ uu (u, σ) = y(u) ⊗ g uu (u, σ). ( 14
)
Finally, the CSS image I c gives in result a multi-scale representation of zero crossing points:

I c = {(u, σ)|κ(u, σ) = 0, u ∈ [0, 1], σ ≥ 0}. ( 15
)
The last stage of the descriptor's construction is the selection of contour maxima locations in CSS image I c . They are used as a representation for a shape.

The similarity between two represented objects is calculated through matching between two sets of maxima, representing them. The method firstly finds any possible changes in orientation in one of the two shapes, then a circular shift is applied in order to compensate their influence. Later the Euclidean distance is calculated between the two sets, indicating the dissimilarity measure.

Discussion on the Experimental Conditions and Obtained Results

The approach described in the previous two sections was experimentally investigated using 60 digital pantomograms, belonging to 30 individuals. For each person two images were applied that were taken at various time points. One image was treated as a template pattern and the second as a test object.

After the pre-processing and segmentation of an image and the extraction of teeth shapes, each obtained contour was represented using the Curvature Scale Space transform. As can be seen on Fig. 3, the CSS shape descriptors obtained for different teeth differ significantly, which confirms the validity of their application in the discussed problem. For the derivation of the similarity between two descriptions the method discussed in the previous section was applied. For further processes the location of each tooth was required, because the matching was performed on the level of particular tooth types -incisors extracted from pantomograms under recognition were matched with incisors from all base images, canines were compared with canines, and so on. The best recognition results were achieved when using the molars. Nevertheless, other tooth types should also be taken into account, since for certain people this type of tooth was not always present.

The indication process was performed separately for each of the four quarters of a pantomogram image. For each of them, the comparison started from the back of the jaw and continued towards the center of the image. If a tooth was present on a particular position in both the test and template sub-images, matching would be performed. The Euclidean distance indicated the template person for a single tooth. The template class with the largest number of indications was selected as the recognized person.

Using the described approach an experiment was performed. As a result of this experiment 20 of the 30 test objects were identified properly, which gives a recognition rate equal to 67%. Although this result may seem far from ideal it is mainly connected with the difficult character of the digital orthopantomogram images.

Conclusions

An approach for automatic people identification using pantomogram images was described in this paper. It consisted of four stages: image enhancement, image segmentation, tooth localisation, and recognition by means of tooth shapes for a particular person. The shapes were represented using the Curvature Scale Space transform and matched with the shapes stored in the template base. The largest number of correctly indicated teeth decided about the selection of the template class that was most similar to the test image.

During the experiment 60 pantomogram images acquired from 30 people were used -half of the images were selected as templates assigned to particular persons, and the other half was used as test objects. This means that two images obtained at various points in time were used for each person. The achieved identification rate was equal to 67%, what is not enough but nonetheless one has to keep in mind the particularly difficult character of the dental orthopantomogram images. Nevertheless, future works on this problem are planned. First, other shape representation techniques can be applied and experimentally investigated, including evaluation of potential hardware implementation [START_REF] Frejlichowski | Analysis of Possible System-Level Hardware Implementation of the Selected Shape Description Algorithms[END_REF]. Second, other methods for the pre-processing of the image can be analysed in order to emphasise the locations of particular teeth before the extraction of their shapes. Finally, other methods for making a final decision about the identified person could be developed and evaluated.

Fig. 1 .

 1 Fig. 1. An exemplary result of the image pre-processing and segmentation, described in this section.

Fig. 2 .

 2 Fig. 2. Examples of the segmented teeth images.

Fig. 3 .

 3 Fig. 3. Examples of teeth contours and obtained CSS descriptors.

Acknowledgements. The authors of this paper wish to thank gratefully to R. Wanat, MSc for his significant help in developing and investigating the described approach.