Duplisea

Louise Savard
email: louise.savard@dfo-mpo.gc.ca

Mathieu Desgagnés
email: mathieu.desgagnes@dfo-mpo.gc.ca

Daniel Duplisea
email: daniel.duplisea@dfo-mpo.gc.ca

SPAM (Sex-structured Pandalus Assessment Model): a stock assessment model for Pandalus stocks

Keywords: Pandalus, stock assessment, age-length structured model, hermaphrodism, Pandalus borealis, Sept-Îles

Introduction

Shrimp fisheries are becoming more and more important throughout the world. According to FAO, 3,120,566 t of shrimp were caught in marine waters in 2008, representing about 60% of the total crustacean production and 3.9% of the world capture fisheries production (FAO 2009). A Pandalidae, Pandalus borealis, is the second most important exploited shrimp species globally, counting for 12.3% of the global shrimp catches (FAO 2010).

Shrimp of the genus Pandalus are exploited in the North Atlantic and the North Pacific oceans mainly by trawl fisheries on large individuals given that shrimp price is often directly linked to its size. This characteristic of the fisheries becomes important because most species of the genus Pandalus (including P. borealis) are protandric hermaphrodites, i.e. they reproduce first as males and then change sex and become females [START_REF] Bergström | The biology of Pandalus[END_REF][START_REF] Berkely | The post-embryonic development of the common pandalids of Bristich Columbia[END_REF]) (Fig. 1).

Sex-change is thus critical for fishery management [START_REF] Fu | The role of sex change, growth and mortality in Pandalus population dynamics and management[END_REF] since females are primarily targeted by the fishery because of their size. Despite a large literature investigating the factors which potentially influence the activation of the sex-change process, there is currently no clear consensus. Density-dependence effects (Koeller et al. 2000a), size or age [START_REF] Koeller | Is size at sex transition an indicator of growth or abundance in pandalid shrimp?[END_REF][START_REF] Wieland | Length at sex transition in northern shrimp (Pandalus borealis) off West Greenland in relation to changes in temperature and stock size[END_REF]), environmental factors [START_REF] Wieland | Length at sex transition in northern shrimp (Pandalus borealis) off West Greenland in relation to changes in temperature and stock size[END_REF]) and evolutionary processes [START_REF] Charnov | Dimensionless invariants for the optimal size (age) of sex change[END_REF] have been studied but none of them alone is a powerful predictor of the yearly variability observed in sex-change. However, sex-change always occurs within a certain length interval, albeit a relatively large one, and consequently, length appears to be a required basis when modelling sex-change. Individuals of P. borealis remain as males for some years. Sex-change occurs in winter and newly transformed females can be identified as primiparous females in the following spring and summer. Primiparous females spawn the following fall. Larvae hatch in spring and reproductive females that survive breeding are then identified as multiparous females.

Multiparous and primiparous females can be distinguished by the presence of sternal spines that start disappearing during summer and are definitively lost little time before their first mating (mating occurs few hours/days before spawning - [START_REF] Mccrary | Sternal spines as a characteristic for differentiating between females of some pandaleidae[END_REF].

Direct determination of age is not possible for P. borealis and age can be approximated by the identification of the first modes in length frequency distributions and by the examination of sex related characteristics [START_REF] Bergström | The biology of Pandalus[END_REF]. Growth of males can be adequately modelled with a von Bertalanffy curve [START_REF] Bergström | The biology of Pandalus[END_REF]; however size at recruitment is variable through time [START_REF] Hvingel | Northern shrimp research in the North Atlantic -state of the art and future research strategy[END_REF] probably as a consequence of varying environmental conditions during the larval and juvenile phases [START_REF] Daoud | Size and temperature-dependent variations in intermolt duration and size increment at molt of Northern Shrimp, Pandalus borealis[END_REF]. Male growth rate decreases for oldest males and consequently the largest male length modes can contain individuals of different ages.

Absence of direct hard pieces such as otolith combined with the difficulty to distinguish cohorts for largest males hinders the development of reliable aging techniques. Tagging techniques is often used but this technique is more difficult to apply to shrimps because of frequent moulting, small size of the individuals, and is generally expensive. Primiparous females are considerably larger than males of the same age because moult frequency increases during the sex change period. Growth of females is then slow and modes of primiparous and multiparous in lengthfrequency samples are often difficult to distinguish.

Despite the economic importance of crustacean fisheries, few analytical methods have been developed for crustacean stock assessments [START_REF] Smith | Methods for stock assessment of crustacean fisheries[END_REF] and more specifically for shrimp stocks. Age structured models have been applied for some crustacean stocks [START_REF] Ices | Report of the Working Group on Nephrops Stocks, ICES CM 2003/ACFM:18[END_REF], including shrimp (ICES 2001), but in most cases they were found inappropriate because of the uncertainties in aging techniques and the difficulties to track the cohorts after sex-change.

Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011

Catch-Survey-Analysis [START_REF] Cadrin | Evaluating two assessment methods for Gulf of Maine northern shrimp based on simulations[END_REF][START_REF] Cadrin | Application of Catch-Survey Models to the Northern Shrimp Fishery in the Gulf of Maine[END_REF] and surplus production models have also been applied [START_REF] Cadrin | A hierarchical approach to determining reference points for Pandalid shrimp[END_REF][START_REF] Nafo/Ices | Report of the NAFO/ICES Pandalus Assessment Group[END_REF] to carry out shrimp assessments. They have the advantage of being simple and rather robust, requiring few data. However, a major limitation is that they do not provide any information on the demographical, and especially sex stock structure which can be considerably altered by commercial harvesting and is essential for fishery management. Consequently, shrimp stock assessments often consist of a descriptive analysis of various indicators and trends (for example commercial catch rates and/or survey abundance indices) through time, sometimes formalized in a traffic light approach (Koeller et al. 2000b).

Length [START_REF] Drouineau | Assessing the impact of discretisation assumptions in a length-structured population growth model[END_REF][START_REF] Drouineau | A length-structured and spatialised model for the Northern stock of European hake (Merluccius merluccius)[END_REF]) and age-length structured [START_REF] Fournier | MULTIFAN-CL: a length-based, agestructured model for fisheries stock-assessment, with application to South Pacific albacore, Thunnus alalunga[END_REF][START_REF] Quinn | A length-based model for hard-to-age invertebrate populations[END_REF][START_REF] Froysa | Fleksibest-an age-length structured fish stock assessment model[END_REF] models have recently raised more and more interest, especially for species such as shrimp, for which growth is poorly known and aging is difficult. Such models would have several advantages for Pandalus stocks. Indeed, this kind of models can be fitted directly to length-structured data, obviating the need of an uncertain and often expensive length to age data conversion. Moreover, a sex and age-length structured model would enable one to model sex-change in a length interval and to provide valuable information on the demographic structure of the stock. A simulation length-based model with explicit sexchange was first proposed by [START_REF] Fu | The role of sex change, growth and mortality in Pandalus population dynamics and management[END_REF] to quantify the importance of some biological processes, especially growth, sex-change and natural mortality, but the model was dedicated to understanding the sex-change process rather than as an assessment method.

Material and methods

Population model

The model has a seasonal time step (Fig. 1) which is required to properly describe male growth, seasonal biological processes (sex-change, reproduction) and fishing patterns (variability in catchability linked to seasonal migrations). stands respectively for the number of primiparous and multiparous females at time step t.

Assumptions regarding length-at-age

Male growth is assumed to follow a von-Bertalanffy growth curve. Consequently, the mean size of the year class y at time step t is:

(1) ( ) ( ) ( )

( ) ( ) t K r y µ L t µ = t µ y y ∆ ⋅ - ∞ - ⋅ - + + exp 1 1 , with r y
µ , the size at recruitment (size at age 1 at the beginning of spring) of the year class and ∆t the duration of time step t.

∞

L is assumed to be known and constant and could be approximated by the maximum length observed for males since growth at age 3 is slow and nearly null at age 4. Following [START_REF] Fournier | MULTIFAN-CL: a length-based, agestructured model for fisheries stock-assessment, with application to South Pacific albacore, Thunnus alalunga[END_REF] and [START_REF] Maunder | A-SCALA: un analisis estadistico de captura a talla estructurado por edad para la evaluacion de las poblaciones de atunes en el oceano pacifico oriental[END_REF], length distribution of an age group at time step t, is assumed to follow a normal distribution with a mean ( ) t µ y and a constant coefficient of variation

Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 cv.

Assumptions regarding sex-change

Quantities related to the sex-change are subscripted with sex. Sex-change is assumed to be length dependent. We assume that the sex-change process is completed at the end of winter (2)

( ) ( ) ( ) ( ) ( ) 2 0, ~ with 1 sex sex sex ε sex 50 sex 50 σ N y ε y e y L = y L ⋅ -
The proportion p(y,l) of males of size l that change sex is given by:

(3) ( ) ( ) ( ) ( )         - ⋅ - y L l R + = l y, p 50sex sex 3 log 2 exp 1 1
However, all males of last age A are forced to change sex.

It is necessary to make an assumption about sex-change during final year Y to provide abundance estimates in year Y+1. Since sex-change is assumed to follow a random walk, it is logical to assume that the sex-change ogive in year Y is similar to year Y-1.

Survival equations

Survival Sr is the result of the natural and fishing mortality. Mean natural mortality M is assumed to be constant over age groups and stages, and is equal at 0.5 year -1 . Yearly deviations are allowed to account for inter-annual variability:

Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011

(4) ( ) ( ) ( ) ( ) 2 0, with M M M ε σ N y ε y e M = y M ⋅ M(y)
denoting the natural mortality in year y.

( )

l t, F m , ( ) t F primi and ( ) t F multi
stand for fishing mortality applied respectively on males of length l, primiparous and multiparous females at time t. Their calculation is explained in the next section. Survival is computed with the following equations:

(

         ∈ - - - ∆t multi F + y M multi ∆t primi F + y M primi ∆t m F + y M m t e = t Sr y t t e = t Sr l t, e = l t 5) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

, Sr

Fishing activity and survey models

Fishing mortality and catches

Quantities related to the commercial fishery are subscripted with a f. Male fishing mortality is assumed to be the product of a selectivity s f (l), a fishing effort E(t) and a year effect q(y) which corresponds to a catchability. Selectivity s f (l) is modelled as a traditional sigmoid function, characterised by two parameters

f 50 L and f R (interquartile range): (6) ( ) ( ) ( )         - ⋅ - f 50 f f L l R + = l s 3 log 2 exp 1 1 Male fishing mortality is equal to: (7) ( ) ( ) ( ) ( ) t E l s y q = l a, t, F f m ⋅ ⋅
Concerning females, it is assumed that selectivity is equal to 1 given that the newly transformed females reach the sizes that are usually fully recruited to the fishing gear ( ShumwayAuthor-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 et al. 1985). Therefore, primiparous females fishing mortality is equal to:

(8) ( ) ( ) ( ) t E y q = t F primi ⋅
Females that bear eggs might have, in certain areas, a more specific behaviour than males or primiparous females. In winter and spring, egg bearing females (named multiparous in the model) aggregate inshore, and consequently have a higher catchability [START_REF] Shumway | Synopsis of biological data of the pink shrimp, Pandalus borealis[END_REF]. A targeting factor multi spring T , is thus introduced into the calculation of multiparous females fishing mortality in spring to describe this interaction:

(9) ( ) ( ) ( ) ( ) ( )       ⋅ ⋅ ⋅ year the of rest the spring in , t E y q t E T y q = t F multi spring multi
Year effect is assumed to follow a random walk to avoid unrealistic inter-annual variations:

(10) ( ) ( ) ( ) ( ) ( ) 2 0, ~ with 1 f f f ε f f σ N y ε y e y q = y q ⋅ -
Given the population numbers and the fishing mortality, expected catches are calculated using traditional Baranov equations:

(11) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )              ⋅ - - ⋅ ⋅ ⋅ - - ⋅ ⋅ ⋅ - - ⋅ ⋅ ∑ ∆t t F + M t F + M t F t N = t C ∆t t F + M t F + M t F t N = t C l t, C = t C ∆t l t, F + M l t, F + M l t, F l , t, N = l t, C multi multi multi multi multi primi primi primi primi primi l m m m m m m m exp 1 exp 1 exp 1 .
Harvest rates H are computed as a ratio of catches in year y over the abundance of the stock at the beginning of that year, i.e. at the beginning of spring, the first seasonal step of a year:

Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), [770][771][772][773][774][775][776][777][778][779][780][781][782][783] The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011

(12) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )                ∑ ∑ ∑∑ ∑∑∑ ∈ ∈ ∈ 1 1 1, t N t C = y H t N t C = y H l al, t N l a, t, C =
(13) ( ) ( ) ( )         - ⋅ - s 50 s s L l R + = l s 3 log 2 exp 1 1
Selectivity is assumed to be equal to 1 for the females which are large enough to be totally recruited to the survey trawl. Abundance indices of primiparous and multiparous females estimated by the model are summed into a single female abundance index IA f to be consistent with surveys that are often conducted in fall when the distinction between the two stages of females is difficult or impossible. Given these assumptions, expected abundance indices are given by:

(14) ( ) ( ) ( ) ( ) ( ) ( )    ⋅ ⋅ ⋅ t N + t N q = IA l , t, N q l s = l t, IA multi primi s fe m s s m .
Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), [770][771][772][773][774][775][776][777][778][779][780][781][782][783] The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011

Observation model and likelihood function

The model is fitted by maximising a likelihood function which is derived from an observation model that describes the uncertainties around observed data.

Contribution of length-frequency data

Length-frequency samples are available from both commercial and survey catches. A sample consists of the frequencies of males in each length class. Following [START_REF] Fournier | estimating growth parameters and age composition from multiple length frequency data sets illustrated using data for southern bluefin tuna (Thunnus maccoyii)[END_REF][START_REF] Fournier | MULTIFAN-CL: a length-based, agestructured model for fisheries stock-assessment, with application to South Pacific albacore, Thunnus alalunga[END_REF] and [START_REF] Drouineau | A length-structured and spatialised model for the Northern stock of European hake (Merluccius merluccius)[END_REF], a robust likelihood is used to account for the samples. Basically, it is based on the property that a proportion in a sample follows a normal distribution, and the likelihood is modified to limit the influence of high and low frequencies. The length-likelihood of a sample samp is:

(15) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) t n + L + l t, f l t, f = l t, v + l t, v l f l t, f e + l t, v = θ | t samp logL samp samp l samp obs samp 1 2 0.1 1 with 10 2 1 log log 2 1 6 2 ⋅       ⋅ -               - ⋅ ⋅ - ∑ - -
where n samp (t) is the number of sampled shrimps and

( ) ( ) ( ) l t, C l t, C = l t, f L = i m m ∑ 1 / and ( ) ( ) ( ) l t, IA l t, IA = l t, f L = i m m ∑ 1 /
are the expected frequency of length class l estimated by the model for a sample from the commercial fishery and from the scientific survey respectively, and f obs (l) the observed frequency in the sample.

Author ) are assumed to follow a gamma distribution. The gamma distribution is a distribution function that is increasingly used in population dynamic models because of its great flexibility [START_REF] Froysa | Fleksibest-an age-length structured fish stock assessment model[END_REF][START_REF] Haddon | Modelling and Quantitative methods in fisheries[END_REF]).

The loglikelihood is given by:

(16) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )         -         -         -         t C t C t C t C + t C t C t C t C + t C t C t C t C ν = θ | t
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )         -         -         t IA t IA t IA t IA + t IA t IA t IA t IA ν = θ | t LogLIA

Parameter estimation

A large number of parameters have to be estimated (Tab. 1). The search domain dimensions quickly increase with the number of years and age groups. Therefore, following [START_REF] Punt | The performance of a size-structured stock assessment method in face of spatial heterogeneity in growth[END_REF], estimation of unknown parameters θ is split into two steps, growth parameters

( ) ( ) { } 0 0 2 , , 1 1 t µ , t µ , µ , , µ cv, K, = θ A r Y r L L
being estimated before the remaining parameters θ 2 (where Y denotes the number of years in the data set, ranging from y 0 to y Y-1 ).

Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783

The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011

Estimating growth parameters

The following procedure is used to estimate the loglikelihood LogL 2 (θ 1 ):

1. given θ 1 , mean size at age ( )

t µ a
is computed for each age group and each time step (equation 1).

2. given θ 1 and given the assumptions regarding the length distribution of an age group, theoretical frequency per length class and per age group is computed.

3. The loglikelihood of a length-frequency sample from commercial or survey catches is computed by a modal analysis with the mean size and standard deviation of the modes kept fixed:

(18) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ∑ ∑                                 ⋅ - - ⋅ ⋅ ⋅ L = l a a a A = a samp samp 2 t µ cv t µ l length t µ cv a t, λ l t, f = θ | t P logL

Estimating other parameters

Given θ 1 , other parameters are estimated by maximising the loglikelihood which is a function of 5 components: (i) the length composition of commercial and survey catches (equation 15), (ii) the total commercial catches (equation 16), (iii) the total abundance from survey (equation 17), (iv) assumption regarding the random walk of L50 sex and (v) assumption regarding the random walk of the year effect q f : (20)

Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783

The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) θ LogL + θ LogL + θ LogL + θ | t LogLIA + θ | t LogLC + θ | t samp LogL + θ | t samp LogL = θ | θ logL f M sex t t surv t com t 1 2 ∑ ∑ ∑ ∑ with: (21) ( ) ∑ ⋅ - y sex sex sex ε σ = θ LogL 2 2 2 1 (22) ( ) ∑ ⋅ - y M M M ε σ = θ LogL 2 2 2 1 and: (23) ( ) ∑ ⋅ - y f f f ε σ = θ LogL 2 2 2 1
Following [START_REF] Drouineau | A length-structured and spatialised model for the Northern stock of European hake (Merluccius merluccius)[END_REF], the loglikelihood is maximised using an evolutionary algorithm [START_REF] Schwefel | Evolution and Optimum Seeking[END_REF]) that provides a starting point to a quasi-Newton algorithm provided in the autodif library (http://www.otter-rsch.com/). The evolutionary algorithm is a stochastic algorithm which is relevant to explore highly dimensional objective functions and which does not require a starting point to be specified unlike traditional quasi-Newton algorithms. The Hessian matrix is estimated by a finite difference approximation and inverted to get the variance covariance and correlation matrices.

Case study: Gulf of St. Lawrence Pandalus borealis stock

Description of the fishery

The northern shrimp (Pandalus borealis) fishery began in the Gulf of St. Lawrence (Fig. 2

Input data

Catch at length is available for males by time step (except in winter) and abundance at length is available for males for the second time step (in summer), both from 1990 to 2008.

Cephalothorax lengths are aggregated into 0.5 mm classes (they are measured at the nearest 0.1 mm) and range from 8.0 to 27.0 mm. In the Sept-Îles area, shrimp larvae hatch in early May [START_REF] Ouellet | Spring oceanographic conditions and northern shrimp Pandalus borealis recruitment success in the north-western Gulf of St[END_REF]. In fall the year after, they are 1½ year old and measure between 8.0 and [START_REF] Daoud | Size and temperature-dependent variations in intermolt duration and size increment at molt of Northern Shrimp, Pandalus borealis[END_REF]. Indeed, the authors found that, at 5 °C which is about the bottom temperature of the area, 20.0 mm males are 4 years old.

Female catch and abundance data are available from 1990 to 2008. Females are split into two stages (primiparous and multiparous females) for the commercial catch for spring and summer.

The separation is not done for the fall season and for the survey. In the Sept-Îles area, spawning occurs in early October. Primiparous females spawn 6 months or so after having changed sex.

Egg bearing females migrate and aggregate inshore in winter and spring. This behaviour has an impact on their availability to the fishing fleet. Therefore, a targeting factor was introduced to the model to account for this migratory behaviour.

Specific parameters values

Some specific parameters values were fixed for the Sept-Îles case study according to expert knowledge (Tab. 2). We set ν com =ν surv =100 which corresponds to a CV of 10%. We set 0.05 = σ f which corresponds to a random walk of coefficient of variation to 5%, given that fishing activity is well known and that catchability is not expected to have changed a lot over the period. Although there is little information on natural mortality inter-annual variability, we choose to fix 0.05 = σ M , to have equivalent weights between natural and fishing mortality deviations. In the absence of information over sex-change, we choose to let more flexibility to the sex-change random walk and fix 0.10 = σ sex .

Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011

The model includes 4 male age groups and is fitted over a 19 years data series. Consequently 111 parameters should be estimated (24 in phase 1 and 87 in phase 2 -Table 1).

Results

Fitting observed data

The model properly fits length-compositions of the male component with modes that can be easily detected (Fig. 3a, Fig. 3b and Fig. 4), demonstrating the adequacy of the growth model and of the use of a time step.

The seasonal pattern of landings was captured as well in the fitting (Fig. 5). Each year, a peak of catches is observed in spring that corresponds to the fishery targeting multiparous females which are aggregated on the hatching grounds. The model poorly fits summer catches for the last two years, probably because the estimates rely on only a few data points. More generally, the fits are a bit poorer for the summer season than for the two other fishing seasons since the late 1990's, perhaps indicating a change in fishing behaviour over the period. Globally, catches have increased through time over the period (Fig. 5) as a consequence of an increase in TAC [START_REF] Dfo | Assessment of shrimp stocks in the Estuary and Gulf of St. Lawrence in[END_REF].

The global trend in total survey abundance indices estimated by the model is consistent with the observations, except for 2003 (Fig. 6). However, a strong year effect in the survey was detected for many species for that year [START_REF] Dfo | Assessment of shrimp stocks in the Estuary and Gulf of St. Lawrence in[END_REF]. Globally, the model tends to smooth the survey signal, which is not surprising given the random walk employed in the model.

Population numbers

Recruitment at age one (Figs. 7 and8) estimated by the model is variable with two periods of

Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 higher recruitments in early 90's and then in early 2000's. Those strong year classes are especially prominent at 15 mm (2 years old) in the length compositions of the 1992, 1999 and 2001 commercial catches (Fig. 3) and at 10 mm (1 year old) in the 1991, 1998 and 2000 survey length compositions (Fig. 4). The abundance of females (primiparous and multiparous) also

shows the same trend with however a few years lag (Fig. 9a). The abundance of females has been gradually decreasing for the last four to five years probably due to the decreasing abundance of males estimated from 2003 to 2007. Recruitment seems to be slightly decreasing in latter years, therefore, a decrease in total abundance is expected in the future (the number of males increased in recent year, but it is probably because of a delayed sex-change).

No obvious stock-recruitment relationship is observed (Fig. 8). The strong year classes are not explained by a high abundance of females and are more likely due to favourable environmental conditions improving larval survival [START_REF] Ouellet | Ocean surface characteristics influence recruitment variability of populations of northern shrimp (Pandalus borealis) in the Northwest Atlantic[END_REF].

Theoretically, the model provides absolute population number estimates, however a rather high correlation is observed between catchability (of both survey and commercial fishery) and recruitment in first year estimates. Consequently, the population number estimates (and therefore the harvest rates) are more likely to be relative abundance estimates.

Natural and fishing mortality rates

The model provides natural and fishing mortality estimates for males and females at each time step (Fig. 9b -male and multiparous female fishing mortality rates are not represented, since they are equal to the primiparous mortality rate multiplied by a constant through time). A high level of both natural and fishing mortality rates at the beginning of the 1990s explains the low abundance of multiparous in the same period. A sudden decrease in natural mortality first,

Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 followed by a more limited decrease in fishing mortality with the combination of high recruitment resulted in a period of relatively high female abundance between 1996 and 2001.

Similarly, the second peak of female abundance in the 2000's happened after a period of high recruitment, lower fishing mortality, and decreasing natural mortality which continued to decrease till 2006. The recent decrease in the abundance of females is probably due to both an increase in the natural mortality rate estimated for the last two years, and a substantial increase in fishing mortality since 2003.

Year effects do not exhibit strong variations through the period; however a model with a constant year effect was rejected by the Akaike Information Criterion [START_REF] Akaike | Information theory as an extension of the maximum likelihood principle[END_REF].

Sex-change

The model relies on a time-varying sex-change-at-length relationship. A time constant function was tested but rejected by the Akaike Information Criterion, demonstrating that interannual variability in sex-change is significant.

The two extreme ogives of sex-change-at-length show that in some years, a significant proportion of males may anticipate sex-change and become females before reaching age 4 or even age 3 (Fig. 10). In other years, exclusively males of age 4 change sex. This may explain why the mode corresponding to age 4 is often hardly distinguished from the mode corresponding to age 3 in length-composition samples. It might also explain why the prediction of the abundance of primiparous females from the abundance of males is uncertain though it would be very important for management. Globally, length-at-sex-change has shown great variations in primiparous length (Fig. 10) which have likely induced variations in female length [START_REF] Dfo | Assessment of shrimp stocks in the Estuary and Gulf of St. Lawrence in[END_REF].

Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4) effectively managed to prevent stock depletion. Moreover, the results of the model on the yearly variations in size at sex-change will probably help to investigate the process and find relevant explanatory covariables. More specifically, density-dependant influence on the variations in size at sex-change may be investigated in the future, since we observed that periods of high or increasing levels in males abundance correspond to periods with low size at sex-change.

Estimating natural mortality yearly deviations was possible because of the absence of fishing in winter due to ice cover. Natural mortality is then not confounded with fishing mortality during this season. The correlation matrix analysis shows that though year to year natural mortality deviations were slightly correlated, they were not correlated to catchability deviations (year effects), confirming that it was possible to at least partially estimate those deviations. The assumption of yearly deviations around a mean seemed sufficient to catch the main trends in the population numbers. However, shrimp natural mortality is known to be correlated to the abundance of important predators such as cod and redfish, which have decreased substantially in many of the Northwest Atlantic shrimp fishing areas. Consequently, it may prove necessary in the future to describe the predation process in the model to clearly distinguish F and M throughout the time series. Moreover, we assumed that M was constant over age and length in the absence of more precise information, more precise mortality-at-length estimates would be a valuable improvement to the model. The present model was specifically developed to assess Pandalus shrimp stocks. Traditional models such as stage-structured models (CSA), surplus production models or age-structured models proved to be poorly adapted to such species because of the absence of aging techniques and the non integration of the sex-change process. In SPAM, the age, length and sex structure provides relevant solutions to both problems. However, a rather large amount of data is required to fit the model, which can be justified by the economic importance of the fishery.

On the whole, results are consistent with the traditional assessment (DFO 2009) confirming a present high abundance level compared to the early 1990's. However, short or medium term projections are possible with SPAM, and consequently, managers may detect more quickly variations in abundance and adapt management measure (TAC) accordingly. More specifically, the detailed male length-structured abundance estimates and sex change at length ogive are valuable information to managers to anticipate more precisely the number of females in the years to come, which are especially targeted by the commercial fishery.

Finally, the model may be used as an operating model in a management strategy evaluation (MSE) framework [START_REF] Butterworth | Experiences in the evaluation and implementation of management procedures[END_REF][START_REF] Punt | Developing management procedures that are robust to uncertainty: lessons from the International Whaling Commission[END_REF][START_REF] Sainsbury | Design of operational management strategies for achieving fishery ecosystem objectives[END_REF].

MSE aims at

  before the beginning of the new biological year. The proportion of males that change sex in a given year y is modelled by a sigmoid function of length, characterised by two parameters follow a lognormal random walk to account for inter-annual variability in sex-change:

  abundance indices Quantities related to the survey are subscripted with an s. Survey abundance indices are assumed to be a relative measure of the stock abundance. Abundance indices are the product of a selectivity s s (l), a catchability q s and numbers of shrimp. We use a sigmoid function, characterised by two parameters s 50 L and s R (interquartile range), to model selectivity for males. Therefore, s s (l) is given by

ν

  com and ν surv equal to the inverse of the square root of the coefficient of variation of the distributions (McCullagh and Nelder 1989).

  λ samp (t,a) estimated for each sample by an Expectation Maximisation algorithm. 4. Repeat steps 1-3 to θ 1 that maximises (19)

Fig. 2 .

 2 Fig. 2. Sept-Îles stock location in the Gulf of Saint-Lawrence.

Fig. 3a .

 3a Fig. 3a. Observed (normal lines) and estimated (bold lines) length-composition of male

Fig. 3b .

 3b Fig. 3b. Observed (normal lines) and estimated (bold lines) length-composition of male

Fig. 4 .

 4 Fig. 4. Observed (normal lines) and estimated (bold lines) male length-composition of survey

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Observed (bold solid line and squares) and estimated (solid line and circles) total

Fig. 8 .

 8 Fig. 8. Estimated stock recruitment relationships: females in autumn against males of age 1

Fig. 9 .

 9 Fig. 9. a) Estimated population numbers obtained from the assessment model : numbers for

  Author-produced version of the article published in Canadian Journal ofFisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 function and (iv) a stage structure (primiparous/multiparous) submodel for females. Contrary to Fu et al. (2001), a simple two stage structure is used for females whereas the male component relies on a length and age structure. The model is then applied to the exploited Pandalus borealis stock of the Gulf of Saint-Lawrence as an illustration showing the relevance of SPAM to assess a commercially harvested shrimp population.

	This paper presents SPAM, a model dedicated to protandric hermaphrodite Pandalids stock
	assessment, that incorporates (i) a length based submodel for male, (ii) a constant growth model
	but time varying size at recruitment, (iii) a time varying length at sex-change relationship

  -produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011

	Contribution of total catches and total abundance indices
	Total commercial catches per time step ( ( ) , t C m	C	primi	( ) , t	C	multi	( ) t	) and total survey abundance
	indices per time step (	IA m	( ) ( ) t IA , l t, fe

  ) in the 1960s. The exploitation is conducted by trawlers in four shrimp fishing areas (SFA) (DFO 2009). In 2008, the Gulf landings reached about 36,000 tons. The Sept-Îles area (SFA 10) is the Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 most productive of the four stocks representing about 41% of the total Gulf catch. The fishery is managed by TAC and the number of fishing licenses is regulated. Fishing season starts on April 1st and closes when the TAC is reached or on the 31st of December at the latest. There is no fishing in winter mainly because the ice cover prevents access to fishing grounds.

	The fishery mainly targets large individuals which are more economically valuable; therefore
	fishermen allocate a large part of their effort in spring on reproductive females on hatching
	grounds. On average, females represented 68% of the Sept-Îles fishing area landings between
	1990 and 2008.

Fishers are required to fill out a logbook, indicating the number of hours fished and a dockside monitoring program ensures control over landings. Both log-book and dockside programs have provided very reliable estimates of fishing effort and catches since 1990. The commercial catches are sampled regularly during the fishing season and a bottom trawl survey is conducted each year at the end of summer. The catch sampling program has been running since 1982 and the research survey has been conducted since 1990. The total number of shrimp measured each year varies from 8,000 to 18,000 for the commercial sampling as well as for the survey.

  Author-produced version of the article published in Canadian Journal ofFisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 12.0 mm. They are between 12.0 and 16.0 mm at 2½ y.o., 16.0 and 19.0 mm at 3½ y.o. and 19.0 and 23.0 at 4½ y.o. The maximum size is 28.0 mm but few males are found at sizes larger than 24 mm. Sex-change occurs after mating (which occurs few hours/day before female spawning), in winter, before they reach the age of 5. Controlled growth experiment studies confirmed this growth pattern

  This study presents a sex, age and length structured model dedicated to shrimp stock assessment. To our knowledge, it is the first assessment model designed specifically for Pandalid shrimps, a species that changes sex and for whom ageing after the sex-change is difficult, if not impossible. Our model is a sex-structured model with a length-based submodel for males and a simple stage-structured submodel for females. To characterise a Pandalid shrimp fishery, modelling the sex-change process from males to females is critical and this length-stage-sexstructured model should thus be considered as at least a way to take this biological process into account. The age-length structure of the male component is required to properly model sexchange given that the size at sex-change varies over time as indicated by the variations in the size of newly transformed females. A more simple stage-structured submodel is used for females because (i) females growth is slow and it is difficult to distinguish modes in length-frequency samples while the presence/absence of sternal spines provide a more reliable discrimination method between primiparous and multiparous females and (ii) females being fully recruited to the fishing gear, their size compositions are not essential for a stock assessment model. The age structure of females can be approximated by the two stages, primiparous (newly transformed females) and multiparous females which can be distinguished by morphological differences.The model demonstrates a high inter-annual variability in length at sex-change. Consequently, the numbers of newly transformed females are difficult to predict from one year to the next. We totally agree with[START_REF] Fu | The role of sex change, growth and mortality in Pandalus population dynamics and management[END_REF] when stating that a better understanding of sex-change would be a valuable improvement for Pandalus stocks management, especially for stocks where females are such a large component of the commercial catches. Females represent the spawning population and constitute the main target of the commercial fisheries and thus they need to be Author-produced version of the article published in Canadian Journal ofFisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011

	, 770-783
	The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4
	doi : 10.1139/f2012-011
	Discussion and perspectives

  The model has been used to assess the Sept-Îles shrimp stock as an illustrative example. The model provides results which are consistent with existing knowledge on the species and the stock Author-produced version of the article published in Canadian Journal ofFisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011this paper was to present an illustrative example. In a formal stock assessment context, those weights should be discussed with experts and adapted to each stock, and sensitivity analysis should be carried out. We think that the model is generic enough to be applied to differentPandalus stocks, if reliable data from both scientific survey and commercial fishery are available. Moreover, the model needs reliable fishing effort data to estimate the seasonal pattern in fishing activity.

	(DFO 2009).
	Weights on sex-change (equation 2 -sex σ ), natural mortality (equation 4 -M σ ) and year
	effect random walks (equation 10 -

f σ ) were fixed rather arbitrarily. However our objective in

Table 1 .

 1 assessing the robustness of management options regarding various sources of uncertainty. It relies on an operating model that is able to simulate realistic population dynamics Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 under various scenarios. The model seems to be able to reproduce Pandalid shrimp stock trends and is flexible enough (especially if one adds time varying natural mortality and size at sexchange) to simulate various plausible scenarios of stock evolution. Unknown parameters and their significance (Y stands for the number of years in the 558 Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011

	559	dataset).
		Type of parameters Parameters							Number of
											parameters
		Growth	K								1
			cv								1
			Mean size at recruitment	y µ ,	r	Y
			Mean size of male age groups (but recruitment)	A-1
			in first time step	µ a	( 0 t	)
		Sex-change	L 50sex (y 0 ) R sex	ε sex	( ) 1 y	ε sex	( ) 2 y	...	ε	( ) 2 -Y sex y	Y
		Natural mortality	ε M	( ) y							Y
		Initial numbers	First year N m (1,a,.) N primi (1) N multi (1)	A+2
			Yearly recruitment N m (y,1,.)	Y-1
		Commercial	Selectivity L 50f , R f					2
		fishery
			Catchability q(y 0 ) ( ) 1 y ε f	ε f	( ) 2 y	... ( ) 1 -Y sf y ε	Y
			Multiparous targetting	spring T	,	multi	for spring	1
		Survey	Selectivity L 50s R s					2
			Catchability q s						1

Table 2 .

 2 Fixed parameters' values for the Sept-Îles case study 560Author-produced version of the article published in Canadian Journal ofFisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 Fig. 1. Illustration of the life cycle of a protandric hermaphrodite Pandalus species and model time-steps (winter, spring, summer autumn) Arrows represent (age or stage or sex) transitions occurring instantaneously between two successive time-steps.

	Parameters Value
	A	4
	L		38
	y 0	1990
	Y	19
	L∞	28
	ν com	100
	ν surv	100
	σ	f	0.05
	σ	M	0.05
	σ	sex	0.1

Acknowledgments

We would like to thank the numerous biologists and technicians who have collected the fishery and survey data, and the PCSH program who funded part of this work. We also like to thank the three anonymous referees for they helpful comments.

Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), [770][771][772][773][774][775][776][777][778][779][780][781][782][783] The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011 Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011

Author-produced version of the article published in Canadian Journal of Fisheries and Aquatic Sciences, 2012, 69(4), 770-783 The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011

the fishing season and b) estimated natural mortality (dotted line) and spring fishing mortality of primiparous females (multiparous mortality is equal to the primiparous mortality multiplied by the targeting factor -solid line). The original publication is available at http://www.nrcresearchpress.com/toc/cjfas/69/4 doi : 10.1139/f2012-011


Rodolphe Vaillant

Loïc Barthe

Gaël Guennebaud

Marie-Paule Cani

Damien Rohmer

Brian Wyvill

Olivier Gourmel

Mathias Paulin

Implicit Skinning: Real-Time Skin Deformation with Contact Modeling

Keywords: CR Categories: I, 3, 7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Animation; Skinning, Mesh deformation with contact

a) (b) Linear blending (c) Our with contact (d) Our with bulge (e) Real finger Figure 1: (a) Real-time animation of the index finger of a hand model composed of 31750 vertices. Two poses are shown in each column: (b) standard smooth skinning with linear blending at 95 frames per second (fps), (c) our method which compensates the loss of volume on top of joints and models contact in the fold at 37 fps, (d) our method with an additional bulge mimicking tissues inflation at 35 fps, and (e) a picture of a real bent finger.

Introduction

One of the main challenges when animating a virtual character is the production of plausible skin deformations at joints (Figure 1). This animation step is called skinning. Popular methods include linear blending skinning (LBS, also called smooth skinning) and dual quaternions [START_REF] Magnenat-Thalmann | Joint-dependent local deformations for hand animation and object grasping[END_REF][START_REF] Kavan | Geometric skinning with approximate dual quaternion blending[END_REF]]. Such geometric skinning techniques give the artist real-time interaction with the 3D model while only requiring limited user input, especially when automatic rigging techniques are used [Baran and Popović 2007;Jacobson et al. 2011;[START_REF] Bharaj | Automatically rigging multi-component characters[END_REF]]. However, the resulting deformation does not maintain a constant volume, and may cause local self-penetrations, thus failing to produce convincing organic-like deformations near joints (Figure 1(b)).

Example-based methods [START_REF] Mohr | Building efficient, accurate character skins from examples[END_REF] and shape interpolation schemes [START_REF] Lewis | Pose space deformation: a unified approach to shape interpolation and skeleton-driven deformation[END_REF] were proposed to achieve more realism while still enabling real-time animation. However, they increase the amount of necessary input and often require tedious work from the artists. Furthermore, they are limited to a given range of deformations, which usually does not include deformations due to contact with other skin parts.

Lastly, physical simulation can also produce appealing deformations [START_REF]Anatomically-based models for physical and geometric reconstruction of humans and other animals[END_REF][START_REF] Teran | Robust quasistatic finite elements and flesh simulation[END_REF]]. However, the amount of computation makes it only suitable for off-line rendering. In addition, muscles and rigid bones of realistic shapes need to be predefined from medical data or by the user. Intermediate methods such as elasticity skinning [START_REF] Mcadams | Efficient elasticity for character skinning with contact and collisions[END_REF]] enables the computation of skin squash and stretch by solving the underlying physical equations. While robust and visually plausible, it still requires several seconds per animation frame.

Our goal is to provide the artist with a real-time skinning technique that automatically produces self-penetration free deformations when skin folds, while preserving the aspect of a rigid bone near joints (Figure 1(c)). We also want to be able to generate subtle effects such as bulges (Figure 1(d)). To this end, we start from an initial mesh deformation produced by standard geometric skinning. The key idea is to approximate the individual mesh parts attached to each skeleton bone using implicit surfaces, and to exploit their ability to be combined using operators such as union, gradient-based blending or gradient-based bulge-incontact [START_REF] Gourmel | A gradient-based implicit blend[END_REF]. Combining these operators enables us to generate plausible skin deformations at joints during motion, which is done by projecting each mesh vertex onto the adequate isosurface, from its smooth skinning position. In addition to a novel skinning approach, our main technical contributions are:

• A specific reconstruction method for mesh parts, based on Hermite Radial Basis Functions (HRBF) [START_REF] Wendland | Poisson disk point sets by hierarchical dart throwing[END_REF]].

• A dedicated blending operator for modeling skinning effects.

• A fast marching algorithm for mesh vertices in a scalar field, used in conjunction with a specific mesh relaxation method avoiding the introduction of large mesh distortions.

The resulting technique generates visually plausible skin deformations in real-time. Our method automatically generates contact surfaces between skin parts, without requiring any collision detection step. Finally, it lies on top of standard geometric skinning and can therefore be easily implemented within a standard animation pipeline.

Related work

Skeleton driven skin deformation was first introduced by the study of a hand model [START_REF] Magnenat-Thalmann | Joint-dependent local deformations for hand animation and object grasping[END_REF]]. This linear blend skinning (LBS) technique combines joint transformations based on skinning weights associated to each bone. Its limitations have been extensively studied, and include the candy wrapper effect (mesh shrinkage when bones twist) and the elbow collapse effect (loss of volume at a bone join).

Many example based techniques were developed as an alternative [START_REF] Lewis | Pose space deformation: a unified approach to shape interpolation and skeleton-driven deformation[END_REF][START_REF] Sloan | Shape by example[END_REF][START_REF] Kry | Eigenskin: real time large deformation character skinning in hardware[END_REF][START_REF] Wang | Multi-weight enveloping: least-squares approximation techniques for skin animation[END_REF][START_REF] Mohr | Building efficient, accurate character skins from examples[END_REF][START_REF] Wang | Real-time enveloping with rotational regression[END_REF][START_REF] Weber | Context-aware skeletal shape deformation[END_REF]]. They achieve realistic results but require the production of several static meshes with the necessary key poses as input in order to provide the desired deformation.

More advanced geometrical skinning methods were introduced to limit the artifacts while keeping their simplicity of use. They relied on matrices interpolation [START_REF] Magnenat-Thalmann | Modeling of bodies and clothes for virtual environments[END_REF]Kavan and Žára 2005], dual quaternion deformation [START_REF] Kavan | Geometric skinning with approximate dual quaternion blending[END_REF], or fast optimization processes [START_REF] Kavan | Automatic linearization of nonlinear skinning[END_REF][START_REF] Shi | Mesh puppetry: Cascading optimization of mesh deformation with inverse kinematics[END_REF]], while possibly allowing bone scaling [Jacobson and Sorkine 2011]. These methods remain however prone to the elbow collapse effect and only produce smooth deformations such as those illustrated in Figures 1(b) and 10. Displacement textures encoded along spline curves were used to generate organic-like deformation [Forstmann et al. 2007]. However, being pre-stored as a texture, the deformation does not adapt to the skinning deformation and does not avoid self-intersections. [START_REF] Kavan | Elasticity-inspired deformers for character articulation[END_REF] proposed to pre-compute optimized skinning weights for linear skinning and dual quaternions at joints to approximate the skin deformations produced by nonlinear variational deformation methods [START_REF] Mcadams | Efficient elasticity for character skinning with contact and collisions[END_REF]. Even though this expensive pre-computation produces more convincing real-time deformations, the solution remains too smooth and does not faithfully reproduce bone joints behaviors where the skin stretches around the joint and where it is crushed due to contact.

Volume preserving skinning methods allow users to correct for volume changes through the generation of extra bulges and/or wrinkles [START_REF] Funck | Volumepreserving mesh skinning[END_REF][START_REF] Rohmer | Exact volume preserving skinning with shape control[END_REF]]. However, this approach is computationally expensive, and contacts between different parts of the skin are not handled. Self-collision can be avoided using a volume preserving algorithm based on vector field integra-tion [START_REF] Angelidis | Kinodynamic skinning using volume-preserving deformations[END_REF]. However this method does not fit into the standard animation pipeline, since it requires some temporal integration. Moreover, the parameters controlling the appearance of the final shape may not be intuitive to set-up.

The idea of using implicit representations in conjunction with meshes for character animation is not new. Metaballs [START_REF] Shen | Interactive shape design using metaballs and splines[END_REF], polygon-based implicit primitives [START_REF] Singh | Implicit surface based deformations of polyhedral objects[END_REF], convolution surfaces [START_REF] Van Overveld | Using the implicit surface paradigm for smooth animation of triangle meshes[END_REF] and ellipsoidal implicit primitives [START_REF] Leclercq | Mixing triangle meshes and implicit surfaces in character animation[END_REF] were used to animate skinned characters, with the ability in the last case to add organic details (muscles) to a coarse mesh representing the characters body. Closer to our work, implicit modeling has been used to improve skinning deformation [START_REF] Bloomenthal | Medial-based vertex deformation[END_REF]]. The shape is encoded through its medial axis, itself animated through skinning, which limits the candy wrapper and collapsing elbow artifacts. However organic effects such as bulges and contact surfaces are not handled by this approach.

Contrary to previous methods, we do not use implicit surfaces to model the shape of the character, but rather to preserve the original mesh properties, produce controllable deformations through compositing operators, and avoid self-collision during deformation. Meanwhile, surface details always remain encoded by the mesh.

Overview

The different steps of our approach are illustrated in Figure 2. As for standard geometric skinning [START_REF] Magnenat-Thalmann | Joint-dependent local deformations for hand animation and object grasping[END_REF][START_REF] Kavan | Geometric skinning with approximate dual quaternion blending[END_REF], we start from a mesh equipped with an animation skeleton defined as a hierarchy of bones, and the associated skinning weights (Figure 2(a,b)). In addition, the mesh is already partitioned with respect to skeleton bones (Figure 2(c)). This input can either be provided by artists or automatically generated. In this paper we use weights computed with the heat diffusion technique of Baran and Popovic [2007].

From these initial settings, we use Hermite Radial Basis Functions (HRBF) [START_REF] Wendland | Poisson disk point sets by hierarchical dart throwing[END_REF]Macêdo et al. 2011] to approximate each part of the mesh with the 0.5-isosurface of a smooth scalar field fi : R 3 → R (Figure 2(d), Section 4.1). After this step, each vertex v of the mesh stores its current field value fi(v) that encodes the details of the mesh. Then, a single field function f is defined from the combination of the fi using either the union [START_REF] Ricci | Constructive Geometry for Computer Graphics[END_REF]], gradient-controlled blending or gradient-controlled bulge operators [START_REF] Gourmel | A gradient-based implicit blend[END_REF], depending in the desired result (Figures 2(e,f), Section 4.2). The key idea is to animate simultaneously the mesh with geometric skinning and the field function f by applying rigid transformations to the underlying field functions fi. The mesh vertices then march following the gradient of f (Figure 2(g)) until they reach their original field value or cross a gradient discontinuity representing a contact surface (Section 5), thus producing the final deformation (Figure 2(h)).

Field function construction

In this section, we explain how the field function that closely approximates the mesh and mimics the skin deformation is computed. First, Section 4.1 presents the individual reconstruction of parts of the mesh associated with a given skeleton using an adaptation of HRBF methods. The construction of the final field function by composition of these HRBFs is then detailed in Section 4.2.

HRBF Reconstruction

Problem setting. Given a bone i and its associated sub-mesh Mi, our goal is to reconstruct a smooth scalar field fi that tightly approximates Mi. In contrast to the standard surface reconstruction problem, our fi needs to satisfy several specific constraints.

Firstly, high field smoothness is essential to avoid unexpected behaviors of the gradient-controlled operators (Section 4.2) and to stabilize gradient-descent projection (Section 5) during animation. Smooth fields also require fewer parameters, which reduces their memory footprint and speeds up the fitting and evaluation steps. In order to reach this smoothness while preventing loss of details of the input shape, our insight is to embed the mesh in the resulting scalar field, rather than having it exactly coincide with the 0.5-isosurface: each vertex vj stores its local field value f (vj) at rest, enabling it to be projected onto its own isovalue during deformation.

Secondly, the field function fi should appropriately close the large holes left by the mesh partitioning near joints as depicted by the blue curve in Figure 3(a).

Lastly, to allow for local compositions of the fi, as well as to speedup the evaluation of f , the field functions should be compactly supported. We also require fi ranges in [0, 1], with 0.5 being the reference isovalue, and use the convention: fi(x) > 0.5 if x is inside, fi(x) < 0.5 if x is outside.

Hermite RBF. Given all these constraints, we propose the use of Hermite RBFs [START_REF] Wendland | Poisson disk point sets by hierarchical dart throwing[END_REF]Macêdo et al. 2011] to reconstruct a global signed distance field di (the mesh being approximated by the zero isosurface). This distance field is then re-parameterized to yield the compactly supported field function fi.

The HRBF is a variant of RBF that enables to seamlessly interpolating Hermite data (set of positions and normals) without having to populate the whole space with RBFs, nor relying on offset constraint points which can produce reconstruction artifacts [START_REF] Shen | Interpolating and approximating implicit surfaces from polygon soup[END_REF]]. In addition, it is scale independent and it robustly reconstructs concavities. More precisely, we seek a distance field di of the form:

di(x) = m k=1 λ k φ( x -v k ) + β β β T k ∇φ( x -v k ) , (1)
where the vertices v k are the HRBF centers, the scalars λ k and vectors β β β k are the unknown coefficients, and φ is a smooth function for which we recommend the polyharmonic spline φ(x) = x 3 . Given a set of m points v k with prescribed normals n k , the 4m unknown coefficients are easily found by solving for the system of

4m equations: di(v k ) = 0 and ∇di(v k ) = n k .
Still remains the delicate choice of the HRBF centers. Since we only want to approximate the input sub-mesh Mi while leaving out the details, a natural choice is to sample the mesh surface with a few samples. To robustly handle arbitrary meshes, we employ a Poisson disk sampling strategy [START_REF] Wendland | Poisson disk point sets by hierarchical dart throwing[END_REF]] (Figure 3(a,b)). In all our tests, targeting for 50 samples has always been sufficient.

Even though degree three polyharmonic RBFs naturally fill the large holes left at the bone extremities (Figures 4(a,b)), they have to be closed in such a way that the common extremities of two adjacent fields slide over each other without introducing gaps, or creating outgrowths as it is the case in Figure 4(b). Indeed, this would result in a poor skin deformation at joint, as illustrated in Figure 4(c).

Ideally, the extremities of two fields at a given joint should be filled with spherical components of the same radius and centered at the joint location, which is impossible in general. Following this observation, we propose to adjust the closure of the distance field as in Figure 4(d) by adding one HRBF center at each extremity along the line supporting the bone i and at a distance sj from the respective joint j (Figure 3). Their normal constraints are aligned with the bone and point outward. The distance sj is set as the distance from the joint to the nearest vertex. This allows us to automatically generate a plausible bone joint deformation as shown in Figure 4(e).

Finally, in order to let the HRBFs produce a very smooth surface, samples which are too close to a joint are automatically removed.

To this end we employ culling planes orthogonal to the bone and placed at a distance h of the extremities. Formally, let b 0 i and b 1 i the two extremities of the bone i, all vertices v k that do not satisfy the following criterion are removed:

h < (v k -b 0 i ) T (b 1 i -b 0 i ) b 1 i -b 0 i 2 < 1 -h . ( 2
)
We found h = 0.05 to be an effective value in all our results.

Re-parameterization. Finally, the compactly supported field functions fi we seek are computed using the following remapping: fi(x) = tr(di(x)), tr being defined as:

tr(x) =    1 if x < -r 0 if x > r -3 16 ( x r ) 5 + 5 8 ( x r ) 3 -15 16 ( x r ) + 1 2 otherwise ,
where r is a value which sets the size of fis compact support. The purpose is to convert the infinite support of the HRBF di into a finite support field function fi, while ensuring smoothness at the border of its support. More precisely, tr is defined as to smoothly map the HRBF values of di which are in [-r, r] onto the interval [0, 1], with tr(-r) = 1 and t ′ r (-r) = 0 (inside the shape), tr(0) = 0.5 (on the surface) and tr(r) = t ′ r (r) = 0 (outside). In order to avoid getting constant values fi = 1 and thus null gradients inside the shape, we set r to the distance between the bone and the farthest sampling point used for reconstruction. Then, fi = 1 is only reached on the skeleton and the support size nicely scales with the size of the reconstructed shape.

Composition into a global field function

Once the individual compactly supported field function fi are computed, we combine them to define a global field function f fitting the whole mesh to be skinned. To this end, we propose the following two level composition strategy. First, each pair of functions fi and fj associated to neighboring skeleton bones are combined using a binary composition operator g k : R 2 → R. The result of this composition is a new field function fk = g k (fi, fj). Then, the final field f is defined as the union of all combined pair of functions (fi, fj) using [START_REF] Ricci | Constructive Geometry for Computer Graphics[END_REF] max operator:

f = max k ( fk ) (3)
This composition strategy allows us to individually control the deformation behavior for different types of joints, as well as to compensate for the lack of advanced n-ary composition operators.

The choice of operators g k depends on the desired deformation. For instance, the union operator fk = max(fi, fj) can be used to generate a fold when the joint articulates while avoiding the elbow collapse effect. This produces the results shown in Figures 1(c), 2(h) and 8(b).

Blending. Smooth blending operators [START_REF] Blinn | A generalization of algebraic surface drawing[END_REF][START_REF] Barthe | Controllable binary CSG operators for "soft objects[END_REF][START_REF] Bernhardt | Implicit blending revisited[END_REF]] such as the sum and bulge-in-contact operators [START_REF] Cani | An implicit formulation for precise contact modeling between flexible solids[END_REF]] have been popular in implicit modeling.

They can be useful to model, respectively, smoother surfaces near joints such as at the elbow of a character, and organic bulging effects with larger contact surfaces as when a finger bends. However, the standard version of these operators would create unwanted bulges at joints, and lead to inadequate surface behaviors when a bone rotates (see Figure 5(b) for smooth blending). The recently introduced gradient-based composition operators [START_REF] Gourmel | A gradient-based implicit blend[END_REF]] overcome these limitations. They interpolate between union and blending depending on a parameter θ, function of the local angle α between gradients of the composed field functions fi and fj at any point p of the ambient space. θ is thus defined as θ(α) and a union is performed wherever θ = θu, a smooth blend where θ = θ b , and an intermediate local blend when θ varies in ]θ b , θu[.

In our case the definition of functions θ(α) is made easier by observing that the angle α between field gradients in the folding region is often close to the value of the rotation angle between the bones associated to fi and fj (see Figure 6(a)).

Let us take the example of joints with quite smooth skin shapes for small angles, as the ones at the elbow or at the knee of a character (as Juna's elbow in Figure 5(a,d)). Setting θ(0) = θu (union), we avoid blending when fi and fj are aligned, i.e. when their gradients are collinear or when their bones are aligned (Figure 5(c) 1 st row), which avoids the unwanted bulge. Then, setting θ(π/4) = θ b (blending), we generate a maximum blend between fi and fj at a position that would correspond to, for instance, a slight flexion of an arm at the elbow, or of a leg at the knee (Figure 5(c) 2 nd row).

Observing that at elbow or knee, the composition should be back to a union (with a sharp crease of the skin at the joint) when α = π/2 (Figure 5(c) 3 rd row) and then remain a union up to the maximal flexion, we also set θ(α) = θu for α ∈ [π/2, π] (Figure 5(c) 4 th row). This specific setting, designed for knees and elbows, leads us to the plot of θ(α) given in Figure 7(a). The corresponding operator equation is computed from the piecewise definition given in [START_REF] Gourmel | A gradient-based implicit blend[END_REF]].

Unfortunately, the original gradient-based blending operator [START_REF] Gourmel | A gradient-based implicit blend[END_REF]] cannot directly be used for our purpose. Indeed, while the 0.5-isosurface of the composed field f roughly represents the skin, its gradient will be used at each animation step for projecting mesh vertices back onto their original iso-values (Section 5). Since geometric skinning shrinks the mesh during deformation, getting adequate gradient vectors inside the shape is especially important. Unfortunately, the existing gradient-based blending operator exhibits an inner depression around the joints, as depicted in Figure 6(b). We thus modified the blending operator to produce an inside-field well suited for the march of the vertices during projection. The equation of our new gradient-based blending operator is given in Appendix A and the improvement in the field is illustrated in Figure 6(c). 1 Bulge-in-contact. Lastly, other types of joints, such as finger joints, almost immediately produce contact surfaces between neighboring skin parts, surrounded by muscular bulges (Figure 1(d,e)). This is best modeled using a gradient-based bulge-in-contact operator as proposed in [START_REF] Gourmel | A gradient-based implicit blend[END_REF]]. This operator allows the interpolation between a union and a contact surface surrounded by a bulge, giving an organic look to the resulting shape by localizing the bulge around contacts only. As for the gradient-based blending operator, choosing specific values for θ allows us to tune the bulging behavior when a joint bends: to avoid the unwanted bulge at the joint in the rest pose, we set θ(0) = θu (union); we then slowly inflate the bulge up to its maximum when the bones are orthogonal (α = π/2), and keep it maximal for larger angles, by setting θ(α) = θc (bulge in contact) , ∀α ∈ [π/2, π], as plotted in Figure 7(b).

2 2 1 2 2 1 3 2 3 3 3 1 2 2 1 2 2 1 3 2 3 3 3 (a) (b)

Surface tracking

During animation, the implicit primitives fi associated with each bone are rigidly transformed and combined as explained in the previous section, resulting in a time-dependent global field f . Directly deforming the mesh by tracking the 0.5-isosurface using the field function derivatives [START_REF] Stam | On the velocity of an implicit surface[END_REF] would result in a smooth distorted mesh after several transformations. Our approach is rather to use the whole region where f is non-zero as an embedding volume for the mesh, used to control its deformations from a standard geometric skinning pose. This enables us to have the mesh track the implicit deformations at low cost, while not losing any surface detail.

Figure 8: Section of the Armadillo's knee of Figure 9 showing the inside part of the mesh. The blue area is the interior side of the mesh and the clear gray its outside. The white lines are the mesh boundaries produced by the cut. As we can see in (a), LBS generates self-intersections in the fold that are converted in (b) into a contact region using our projection.

Dual quaternions [START_REF] Kavan | Geometric skinning with approximate dual quaternion blending[END_REF]] are used to compute an initial position of the mesh vertices at each animation step. The challenge then consists in designing a fast projection operator to correct vertex positions while accounting for both the detailed mesh shape at rest, and the current implicit deformation modeled by f . Our solution makes use of projection in gradient field directions, tangential relaxation, and local Laplacian smoothing, as detailed below.

Vertex projection. During deformation, we project each vertex vi back onto its original iso-value isoi (with isoi = f (vi) in the rest pose), enabling us to account for both the detailed mesh shape and for the current deformation. This is done using Newton iterations, which are both faster and more robust than the use of a constant step size:

v i ← v i + σ (f (v i ) -isoi) ∇f (v i ) ∇f (v i ) 2 . ( 4
)
Using σ = 0.35 is an effective compromise between speed of convergence and accuracy.

In order to prevent any self-intersection of the final surface (Figure 8(a)), as well as to preserve the regularity of the mesh, the projection has to be stopped at the contact surface between different skin parts (Figure 8(b)). Contact surfaces correspond to gradient discontinuities of the scalar field. We detect the later using the following conservative heuristic: during projection, we keep track of the angle between the gradient at two consecutive iterations, and stop when it is greater than 55 • , thus considering we are on the contact surface (see Figures 2(h) and 8(b)).

Tangential relaxation. Since the initial vertex positions provided by geometric skinning might be far away from their final position, simply using the previous method could introduce high distortions of mesh faces, and at the extreme cause self-intersections. In order to both minimize them and improve the march of the vertices in the gradient field, we interleave projection steps with tangential relaxation steps, which move each vertex towards the weighted centroid of its neighbours. More precisely, given a vertex vi, let qi,j be its one-ring neighbors projected onto its tangent plane. As a preprocess, we compute the barycentric coordinates Φi,j such that vi = j Φi,j qi,j, using the mean value coordinate technique [START_REF] Hormann | Mean value coordinates for arbitrary planar polygons[END_REF]. Each relaxation step moves the vertices tangentially to the local surface using: where µ ∈ [0, 1] controls the amount of relaxation. While by construction, this has no effect on the rest pose, it improves the distribution of the vertices after deformation. A relaxation step is applied after each pair of vertex projection steps with an adjusted value of µ, such that the displacement of the vertices decreases with the distance to their target isosurface. This is done by setting:

vi ← (1 -µ)vi + µ j Φi,j qi,j , ( 5
) (a) (b) (c) (d)
µ = max 0, 1 -(|f (vi) -isoi| -1) 4 . ( 6
)
Laplacian smoothing. Even though we assume smooth inputs, a deformation performed by the union operator introduces sharp features at the boundaries of contact regions (Figure 9(b)). The latter should be smoothed out in order to mimic realistic skin deformations. We remove these high frequencies by locally applying Laplacian smoothing to the final mesh, at and around contact regions:

vi ← (1 -βi)vi + βi ṽi , ( 7
)
where ṽi is the centroid of the one-ring neighborhood of vi, and βi controls the amount of smoothing. This is only done for mesh vertices marked as in contact regions, so that surface details are preserved elsewhere. To this end, we set βi to 1 for vertices stopped at a gradient discontinuity, and to 0 for the others, and smooth the βi values over the mesh by diffusion, prior to applying Laplacian smoothing (Figure 9(c)). As illustrated in Figure 9(d), this effectively smoothes out the mesh oscillations at sharp features, resulting in organic-looking shapes. Note that this smoothing step has to be performed after the projection step. Therefore it cannot be combined with the relaxation that is interleaved with projection during the march of the mesh vertices.

Implementation and results

Implementation. All our results were generated on a Intel Core i7 3770K at 3.5GHz, with 16GB of memory and a Geforce 680.

Our CUDA implementation makes extensive use of GPU parallelism. In particular dual quaternions, projection, tangential relaxation and localized Laplacian smoothing are all parallelized over Linear blending Table 1: For each model, from left to right: the number of vertices, the number of skeleton bones, the memory consumption for the storage of HRBFs in MB, the frame rates in fps when animating a single joint using our technique, and the frame rates in fps when moving simultaneously most bones in a real animation (poses are shown in Figure 11) with dual quaternions and our technique.

Dual quaternions Implicit

(b) (c) (d) (e) (f) (g)
the vertices. The individual field functions (Section 4.1) and composition operators (Section 4.2) being both compactly supported, are respectively sampled into 3D textures of resolution 32 3 (computed from the HRBF equation in 15 ms) and 128 3 with trilinear interpolation, leading to very fast evaluations of f and g k on the GPU. Our models are composed of 20 to 70 bones, requiring 10Mb to 35Mb of memory for storing the field functions (Table 1).

The performances of our method for animating the different models used in the paper are summarized in Table 1. The frame rates only stand for the deformation time. As we can see, the performance of our technique mainly depends on the number of deformed mesh vertices that are to be projected on the implicit surface. However, the parallel implementation performed on recent GPUs allows us to animate simultaneously most bones of models composed of thousands of vertices at more than 60 fps and models composed of tens of thousands of vertices at more than 12 fps (examples of animation poses are shown in Figures 10 and11). Even highly detailed models composed of more than 170, 000 vertices are still animated at several fps. allows, in addition, the creation of a smooth skin deformation inside the bend region for small rotation angles, before generating the contact for angles greater than π/2. This is especially useful where the generation of a hollow or a contact for small angles is unrealistic. More subtle and also more realistic, for angles larger than π/2, the contact area is kept farther from the joint (Figure 12(b)).

Finally, the bulge-in-contact operator can be used to mimic the tissues/fat bulging when fingers bend. This is illustrated in Figure 1(d) on phalanx deformations during the hand animation. In all these animations, smooth deformations and contact are adequately generated while the mesh details are preserved, as illustrated on the Armadillo in Figure 11, its thigh and calf in Figure 9 and on Juna's shoulder in Figure 12(d).

Our technique also addresses the candy wrapper effect when bones twist. This is illustrated on a twisted cylinder in Figure 13 where we can see that the deformation is adequately corrected with our operators, even when the bones bend. We illustrate our result in Figure 13(c) using gradient-based blending as it is the most challenging operator to handle with our projection method. If not considering the composition operator that is left as a free parameter able to perform a union, a blend or a bulge-in-contact, the deformation is controlled by preset parameters. In our technique, the resulting shape is not very sensitive to a slight variation of these parameter values and modifications are predictable. Thanks to the interactive feedback provided by our method, we could easily set these parameter values experimentally in interactive sessions, as illustrated in the accompanying video with the adjustment of the bulge size for a finger. Also, when reconstructing the mesh with implicit surfaces, we are seeking for a smooth approximation of mesh parts, independently from the part size and geometry. Fifty samples were sufficient for all our models, which are of different resolution and include different levels of detail. Automatically modulating this number according to some given approximation error may however be a useful extension.

In our tests, complex joints such as shoulders may require the user to add or remove a few sample points used as HRBF centers at the vicinity of the joint. Here, the number of points to be added or removed mainly depends on the appropriateness of the input partitioning solution (Figure 2(c)). In our framework, this action is again performed with real-time feedback.

By default, the composition operator is set to the union for all the model bone joints, as was done for the Dana model (Figure 10). The result is already a more plausible deformation than the one provided by LBS, dual-quaternions or the more recent elasticityinspired deformer [START_REF] Kavan | Elasticity-inspired deformers for character articulation[END_REF] thanks to the rigid bone aspect generated outside the fold and to the contact produced inside. The deformation can then benefit from the gradient-based blend and bulge-in-contact implicit composition operators to be improved at, for instance, elbows, shoulders and knees for the former and fingers for the later. These two pre-set gradient-based operators are provided to the user: he only has to select between a sharp (union), a smooth (blending) or a bulging (bulges-in-contact) deformation at each joint.

Discussion and limitations

Influence of the initial geometric skinning solution. At each animation step, the mesh is first deformed using a standard geometric skinning, before being corrected through the projection of its vertices to the adequate isosurfaces. Therefore, the quality of our results is affected by the choice of the initial skinning solution. In particular, our correction gives better results if the diffusion of the bone influences is large enough to produce a smooth initial deformation. The geometric skinning used should also avoid the apparition of deep self-intersections which would send some mesh vertices beyond the bones. In that case, these vertices would be projected on the opposite surface, as illustrated with the union operator in Figure 14. In the first row of the figure, the smoother dual quaternions deformation enables an adequate correction, while the input LBS deformations in the second row generate deep selfintersections, causing a failure of our method. Also, in extreme poses, even though contact is handled, the deformation may be inaccurately corrected: some mesh vertices may cross a contact and be stopped, while they were supposed to be projected on the visible skin part, outside the fold.

Smooth skin parts. The neck of a giraffe, the elephant's trunk or the arms of an octopus smoothly deform when they are animated. This effect is not handled by our method which focuses on more rigid bone-like joints. In this case, the smoothing property of LBS near a joint, usually seen as an artifact, may produce a result closer to the expected shape.

Time independence.

At each animation step, the initial mesh vertex positions are computed from the rest pose after a geometric skinning deformation. Therefore, the final vertex positions at a given animation step do not depend on their positions at the previous step. The consequences are two-fold. On the one hand, the solution might be subject to flickering. In practice, as soon as the mesh resolution is dense enough to capture the local shape at joints, the continuity in time of the initial solution computed with geometric skinning coupled with the smoothness of the field functions prevent this effect. On the other hand, this eases the integration of our approach in standard animation pipelines where frames are computed in parallel.

Mesh resolution.

As mentioned in the previous paragraph, the mesh resolution should be dense enough to capture features such as contact, blend or bulge. When the mesh has too coarse a resolution, very few vertices lie on the different features introduced by our skinning technique at joints. Thus, they are not appropriately represented and the relaxation followed by the local smoothing (at the end of the projection) are likely to degrade the vertex distribution, especially if the mesh is very irregular. This can make the mesh deformation unrealistic and discontinuous in time.

Conclusion and future work

We presented a novel approach for skinning virtual characters, which operates in real-time. Based on the embedding of the skin mesh into a deformable implicit volume, made of parts rigidly attached to each bone, it does not suffer from the loss of volume inherent to standard geometric skinning methods. Choosing the way the implicit parts are combined enables us to tune transitions between smooth and sharper skin shapes when a joint bends, to generate contact surfaces between neighboring skin parts, and to drive organic bulging effects. As a result, our approach achieves visually plausible deformations of the different joints in a character's body, even for extreme bending angles. Since no optimization or collision processing steps are required, our approach is robust and efficient. Lastly, computations being independent from one frame to the next, the method perfectly fits into standard animation pipelines.

With our technique, the union operator is used to generate contacts between non-neighboring reconstructed skin parts, as done in the teaser image. In the future, we would like to investigate more advanced field combination schemes to model a wider set of deformations due to collisions between these non-neighboring parts. This may require extending the binary gradient-based composition operators to n-ary compositions. The field values along contact surfaces could also be used, either to output contact forces to be applied to skeleton bones, or simply to detect that an angular limit is reached at the joint.

We would also like to develop more general combination operators, enabling us to capture both smooth blends for small angle values, and contact surfaces surrounded by bulges for larger bends. The fact that gradient-based operators can be edited in real-time by playing with the function θ should ease this tuning. We could also imagine an expert mode where users could edit and control skin deformation at joints through direct interaction, without requiring any knowledge on the theory of our method: after selecting the deformation (sharp, smooth or bulging) for a joint, they could adjust its parameters by directly articulating the input mesh into a specific bending angle, tuning shape composition there, and looking at the resulting animation in real-time.

Finally, inspired by the work of [START_REF] Rohmer | Exact volume preserving skinning with shape control[END_REF] and [START_REF] Barthe | Extrusion of 1D implicit profiles: Theory and first application[END_REF], implicit skinning could be directly tuned from profile curves sketched by the user, and depicting the desired skin shape. These profiles curves, possibly including wrinkles such as those that appear when a wrist articulates, would drive the generation of specific gradient-based compositions.

Figure 2 :

 2 Figure 2: Overview. (a) An input mesh with its animation skeleton, (b) deformation weights at a joint, and (c) mesh segmentation. (d) Implicit surfaces computed as 0.5-isosurfaces of HRBFs approximating each part of the mesh. (e) Composition with a union operator, and resulting shape. (f) Each vertex is assigned to a field value in the rest pose. (g) During animation, both the mesh and the field functions are transformed and mesh vertices march along the field gradient until they reach their individual iso-value, unless they reach a gradient discontinuity modeling some contact region. (h) This produces the final skinned mesh. In (f), (g) and (h), the isosurfaces of the field function are plotted in a vertical plane centered on the skeleton. In blue, the outside values (f < 0.5) and in red the inside values (f > 0.5).

Figure 3 :Figure 4 :

 34 Figure 3: Reconstruction of a mesh part (a phalanx of the hand shown in Figure 2). (a) HRBF nodes are uniformly spread on the mesh surface, and two additional nodes in red are aligned with the bone to close the holes. (b) The resulting implicit surface, and (c) its associated compactly supported scalar field.

Figure

  Figure Arm of the Juna model (Figure 12(a)) and its partitioning. (b,c) Plot of a plane section of the field function f produced by two different compositions of the two field functions reconstructing the arm around the elbow at different rotation angle: (b) With a standard blending and (c) with our new gradient-based blending. (d) The resulting mesh after vertex projections when the gradientbased composition shown in (c) is used.

Figure 6 :

 6 Figure 6: (a) The blue arrows are gradients of the upper right limb field function and the red arrows are gradients of the lower left limb field function reconstructing the arm of Figure 5(a) after an elbow flexion of π/4 (fields are composed with a union). On the bottom, α = π/4, and on the top, gradients are collinear (α = 0). (b) Illustration of the field distortions introduced by the gradientbased blending operator inside the shape, around the joint. (c) Our gradient-based blending operator avoiding this undesired depression.

Figure 7 :

 7 Figure 7: Plot of our functions for θ(α), (a) to get skinning with some blending (around π/4) at an elbow or a knee, and (b) to mimic the inflation of tissues in contact situation.

Figure 9 :

 9 Figure 9: Animation of the Armadillo's knee (a). (b) The result of the projection without smoothing. (c) In red, the smoothed area where βi > 0 and (d) the result after a few steps of Laplacian smoothing.

Figure 10 :

 10 Figure 10: Dana model in a break-dance pose. (a) From left to right, the model is deformed with LBS, dual quaternions and our technique. Note the very smooth deformation with visible loss of volume produced by the LBS (left). (b) Close-up on Dana's knee deformed with dual quaternions and (c) after correction with our method using the union operator. (d) The difference is shown by superposing the two surfaces in (b) and (c). The dual quaternion deformation is in flat-blue with a dark blue silhouette, and our improved solution is transparent gray with wire-frames and a brown silhouette. The overgrowth produced by dual quaternions is corrected and contact surfaces are generated at knee and crotch. (e, f, g) Same illustrations on Dana's elbow.

Figure 11 :Figure 12 :

 1112 Figure 11: Different models skinned with our technique. Number of bones and vertices, memory and frame rates are given inTable 1.

Figure 13 :

 13 Figure 13: Illustration of our method for modeling smooth skin surfaces at joints (new gradient-based blending operator), when the bones respectively twist (top row), or twist and bend (bottom row). (a) LBS, (b) dual quaternions, (c) our technique with gradientbased blending.

Figure 14 :

 14 Figure 14: field of a bent cylinder, reconstructed with a union operator. In red, the part inside the 0.5-isosurface, and in blue the outside. The mesh is shown in transparency. (a) The mesh deformed using dual quaternions is (b) adequately corrected by our method. However, (c) the mesh deformed with LBS, has sets of vertices crossing each other, and located beyond the bones. (d) It leads to many projections in the wrong direction where our method fails.

Acknowledgments

This work has been partially funded by the IM&M project (ANR-11-JS02-007) and the advanced grant EXPRESSIVE from the Eu-Research council. Partial funding also comes from the Natural Sciences and Engineering Research Council of Canada, the GRAND NCE, Canada and Intel Corps. Finally, this work received partial support from the Royal Society Wolfson Research Merit Award.

A Gradient-based blending operator

In this Appendix, we give the equation of our gradient-based blending operator. It has been specifically designed for improving the march of the mesh vertices in the field during projection. Starting from the original gradient-based blending operator presented in [START_REF] Gourmel | A gradient-based implicit blend[END_REF], we modify the inside part such that the blend decreases between fields values 0.5 and 0.7. Deeper inside, when fi > 0.7 or fj > 0.7, the operator is a max (union) (Figure 6(c)). This leads to the following definition for our new blending operator g b :

in which ḡC is defined as :

and the function k θ (C) is defined as:

In practice, equation 8 does not need to be explicitly solved. We rather pre-compute the values of g b in a 3D texture using the method in [START_REF] Gourmel | A gradient-based implicit blend[END_REF]]. Note that this procedure requires the definition of a function s : R → R. In our case, we simply use s = 1.