Vladimir Fux

Patrick Maillé

Matteo Cesana

Price Competition between Road Side Units Operators in Vehicular Networks

Vehicular networks, besides supporting safetyoriented applications, are nowadays expected to provide effective communication infrastructure also for supporting leisureoriented application including content sharing, gaming and Internet access on the move. This work focuses on Vehicle to Infrastructure (V2I) scenarios, where multiple content providers own a physical infrastructure of Road Side Units (RSUs) which they use to sell contents to moving vehicles. Content provider/RSU owners compete by adapting their pricing strategies with the selfish objective to maximize their own revenues. We study the economics of the price competition between the providers by resorting to game theoretic tools. Namely, we formalize a simultaneous price game among the operators further studying the existence of Nash equilibria and their related quality in terms of Price of Anarchy and Price of Stability. The proposed game model is finally used to assess the impact onto the game equilibra of several practical factors including the vehicles' willingness to pay, the traffic densities, and the configuration of the physical networks of RSUs.

I. INTRODUCTION

The constant increase in the number of cars traveling along the roads worldwide calls for effective means to improve the road safety and the efficiency of the overall transportation infrastructure. To this end, the research community, the industries and the governments all over the world are investing much of their efforts and money on the development of integrated Intelligent Transportation Systems (ITS) based on wireless communication networks allowing vehicles, equipment on the road, service centers and intelligent sensors to exchange information in a prompt and cost effective way. In this scenario, vehicles are geared with wireless communication hardware, often referred to as On Board Units (OBUs), to support communication with other vehicles (Vehicle-to-Vehicle, V2V) and with road infrastructure (Vehicle-to-Infrastructure, V2I). In this last case, the devices composing the roadside infrastructure are often called RoadSide Units (RSUs).

A broad classification of the applications which are enabled by vehicular networks can be found in [START_REF] Hartenstein | A tutorial survey on vehicular ad hoc networks[END_REF] where a distinction is made between applications targeting safety, transport efficiency, and information/entertainment. Safety applications include, as an example, collision warning services, transport efficiency application may include lane merging assistance, and navigation services, whereas information/entertainment application range from file sharing among vehicles to Internet access on the move.

In this work, we focus on the vehicle-to-infrastructure (V2I) communication paradigm for VANETs to support content distribution to moving vehicles. Namely, we consider the case where multiple content providers coexist and compete in a given geographical area. Each content provider owns a physical infrastructure of RSUs which she uses to sell contents to moving vehicles. Content provider/RSU owners compete by adapting their pricing strategies with the selfish objective to maximize their own revenues. In such a scenario, we ask ourselves the following simple question: if competing providers wish to select the pricing strategy in order to provide or collect data to/from passing vehicles, what kind of strategies should they follow? The answer is far from being trivial as it predictably depends on several factors including the vehicles' willingness to pay, the traffic densities, the configuration of the physical networks of RSUs, and the strategic interaction among the content providers.

We tackle this problem by considering a basic model with a duopoly of competing content providers. We study the economics of the competition between the two providers by resorting to game theoretic tools [START_REF] Osborne | A Course in Game Theory[END_REF], [START_REF] Maillé | Telecommunication Network Economics: From Theory to Applications[END_REF]. We take a constructive approach by first analyzing the best response of strategy for one operator when the pricing strategy of the other one is given. Then, we extend the analysis to the case where both operator compete simultaneously to set their pricing strategy. In this last case, we formally study the existence of Nash equilibria for the duopoly pricing game and their related quality in terms of Price of Anarchy.

The manuscript is organized as follows: Section II reviews the related literature in the field further highlighting the major novelties of the present work; in Section III, the reference network scenario and model are described. Section IV provides the analysis of the best response of one operator when the pricing strategy of the other is given, whereas the simultaneous pricing game is analyzed in Section V. Section VI concludes and indicates some directions for future work.

II. RELATED WORK

The design of efficient V2I and V2V networks has already attracted much attention within the research community. Most of the work generally targets the design and optimization of communication protocols to be used in vehicular networks. As an example, the optimization of V2I segment is targeted in [START_REF] Zhang | On scheduling vehicle-roadside data access[END_REF] where the focus in on uplink and downlink packet scheduling techniques. Along the same lines, Yang et al. study in [START_REF] Yang | A multihop peer-communication protocol with fairness guarantee for IEEE 802.16-based vehicular networks[END_REF] the applicability and performance of IEEE 802.16 for the communication between groups of vehicles and an RSU.

V2V communications are addressed in [START_REF] Borgonovo | ADHOC MAC: New MAC architecture for ad hoc networks providing efficient and reliable point-to-point and broadcast services[END_REF], [START_REF] Lee | Dissemination and harvesting of urban data using vehicular sensing platforms[END_REF], [START_REF] Yang | A vehicle-to-vehicle communication protocol for cooperative collision warning[END_REF]. In [START_REF] Borgonovo | ADHOC MAC: New MAC architecture for ad hoc networks providing efficient and reliable point-to-point and broadcast services[END_REF] a Medium Access Control (MAC) protocol is proposed to support reliable communication among vehicles. The work in [START_REF] Yang | A vehicle-to-vehicle communication protocol for cooperative collision warning[END_REF] proposes a protocol framework to support the dissemination of warning messages in V2V, whereas the use of V2V communications to support proactive data monitoring in urban environments in studied in [START_REF] Lee | Dissemination and harvesting of urban data using vehicular sensing platforms[END_REF].

In the field of V2I networks, besides the work on protocol design/optimization, it is worth mentioning the research field targeting the optimal design of the roadside infrastructure. In this case, the goal is to optimize the deployment of the RSUs with respect to specific objectives which are generally related to the coverage ratio of vehicles. Trullols et al. [START_REF] Trullols | Planning roadside infrastructure for information dissemination in intelligent transportation systems[END_REF] propose three formulations for the the deployment problem as a Maximum Coverage Problem (MCP), Knapsack Problem (KP), and Maximum Coverage with Time Threshold Problem (MCTTP), respectively; heuristics based on local-search and greedy approaches are then introduce to get suboptimal solutions. Along the same lines, Cavalcante et al. [START_REF] Cavalcante | Roadside unit deployment for information dissemination in a VANET: An evolutionary approach[END_REF] focus on the Maximum Coverage with Time Threshold Problem (MCTTP) and propose a genetic algorithm to solve it. Yan et al.. [START_REF] Yan | Access points planning in urban area for data dissemination to drivers[END_REF] study the very same RSU deployment problem in case the candidate sites for deployment are limited to the intersections between crossing roads. The interested reader may refer to [START_REF] Barrachina | Road side unit deployment: A density-based approach[END_REF] and references therein for a more comprehensive description on the general problem of RSU deployment. Different from the aforementioned work which assumes one central entity to optimize the RSU deployment, [START_REF] Filippini | Non-cooperative RSU deployment in vehicular networks[END_REF] studies the competitive scenario where different network operators compete in the deployment of their respective RSUs by resorting to a noncooperative game. Spatial positioning games are also proposed in [START_REF] Altman | Spatial SINR games combining base station placement and mobile association[END_REF] for generic wireless access networks.

Game theory has been used to evaluate the strategic interaction between the different agents in vehicular networks [START_REF] Niyato | Game-Theoretic Models for Vehicular Networks[END_REF]. In [START_REF] Niyato | Competitive wireless access for data streaming over vehicle-to-roadside communications[END_REF], the authors introduce a stochastic game among OBUs which compete to get service from shared RSUs. Nyiato et al. propose in [START_REF]A unified framework for optimal wireless access for data streaming over vehicle-to-roadside communications[END_REF] a hierarchical game framework to capture the competition of different actors; besides OBUs and RSUs, the concept of Transit Service Provider (TSP) is used to model an entity which manages groups of vehicles and is in charge of minimizing the total cost to support streaming application to its vehicles while meeting the application QoS requirement. The available bandwidth at each RSU can be split in reserved bandwidth and on-demand bandwidth. OBUs make short-term decisions between on-demand and reserved bandwidth (if available), TSPs decides what kind of bandwidth split to purchase from different RSUs along the road, whereas Network Service Providers (NSPs) owning RSUs set their price for on-demand bandwidth to maximize their revenues. Differently, in [START_REF] Saad | Coalition formation games for distributed cooperation among roadside units in vehicular networks[END_REF] a coalition formation game among RSU is analyzed, with the aim of better exploiting V2V communications for data dissemination.

The matter of pricing in generic wireless access networks is largely debated in the literature. Reference [START_REF] Niyato | Competitive pricing in heterogeneous wireless access networks: Issues and approaches[END_REF] provides a nice overview on pricing problems in wireless networks, and further analyze a specific case where two wireless Internet service providers compete on prices, one owning a WiMAX-based infrastructure and the other running a WiFi-based infrastructure. Differently from previously mentioned literature, in this work we focus on price competition between network operators for V2I networks, which is, to the best of our knowledge, a novel issue. Even if V2I networks bear some similarities with generic wireless access networks, there are distinctive features which make the pricing problem worth analyzing; in generic wireless access networks, the network operator competition is generally on the "common" users, that is, those users which fall in the coverage area of the competing network providers. In other words, there is actually a competition only if the coverage areas of the network providers (partially) overlap as in [START_REF] Niyato | Competitive pricing in heterogeneous wireless access networks: Issues and approaches[END_REF]. Users themselves tend to choose the network operator which maximizes some quality measure as in [START_REF] Fux | Incentivizing efficient load repartition in heterogeneous wireless networks with selfish delay-sensitive users[END_REF]. On the other side, in V2I networks competition may arise due to vehicles mobility even if the coverage areas of competing RSUs are not overlapping, since if a RSU does not serve a moving vehicle in its own coverage range, the very same user can be served later by competing operators.

III. MODEL

A. Usage scenario

We consider two Internet access providers (labeled by 1 and 2), competing to attract users on a stretch of a highway. They offer the possibility to access the Internet through Road Side Units, which allows cheaper or better QoS than the other available cellular networks 1 . We assume that each provider has already deployed one RSU -on different locations along the road-, and that both RSU are identical; we denote their individual goodput (or capacity) by c. (Note that this model easily extends to the case when providers own disjoint "connectivity regions", each one made of several RSUs and with total service capacity c.)

Since both providers' RSU are at different locations, vehicles taking the road in one direction first enter the coverage area of Provider 1's RSU, while those traveling in the opposite direction first see Provider 2. We denote by ρ j , j = 1, 2 the average number of commuters per time unit that first enter Provider j's coverage area; they will cross the competitor's coverage area afterwards (since we are considering only one road).

Each user wants to download data files, for an average volume per user (assumed independent of the travel direction) normalized to 1 without loss of generality; the potential demand (in volume) from users seeing Provider j first thus also equals ρ j . In this paper, we treat those average loads as static values, i.e. we do not model the time variations of the load. Moreover, we assume that the coverage area size of RSUs and the vehicles' speed do not constrain the transfers: if a RSU's capacity exceeds its (average) load, all requests are successfully served. Each provider j = 1, 2 chooses the (flat-rate) price p j to charge for the connection service. To model heterogeneity among users, we assume that only a proportion w(p) of users accept to pay a unit price p for the service (this being independent of the download volume). As a result, if Provider j sets his price to p j , the users who first enter Provider j's service area generate a demand (again, per time unit, and treated as static) of w(p j )ρ j . Note that we are assuming here that users do not try to anticipate the price set by the next provider: when a user first sees an RSU access offer, she responds to it as if there were no other RSU afterwards.

Figure 1 summarizes that scenario in terms of demand flows. The total potential demand (volume per time unit) ρ j from users seeing Provider j can be decomposed into:

1) users accepting the price p j and being served by Provider j; 2) users accepting the price p j and being rejected due to the RSU capacity limit (forming a spillover flow ρ sp j heading to the competitor's RSU); 3) and users refusing the price p j (forming a flow ρ ref j heading to the competitor's RSU). The two latter flows then enter the coverage area of the competing provider, where they can be served or not. In the latter case, we denote the corresponding (unserved) demand by ρ us j . Note that we assume users keep the same willingness-topay for the service when they enter the second RSU coverage area.

B. Mathematical formulation

We now give analytical expressions for the different demand components, using the RSU capacity c and the willingnessto-pay function w(•). In the whole paper, w(•) is assumed continuous and non-increasing, and such that w(0) = 1 and w(p max ) = 0 for some p max > 0. If the quality of the alternative cellular access (say, 4G) is sufficient, the price p max may be interpreted as the unit price for that cellular service: above p max , users have no interest to use an RSUbased access.

The demand submitted to Provider j comes from three different types of users:

1) those seeing Provider j first, and accepting to pay the proposed price p j , hence issuing a total demand w(p j )ρ j ;

2) those seeing Provider k = j (the competing provider) first, who refused to pay p k but would accept the price p j , forming a total demand level (smaller than ρ ref k , and null when p k ≤ p j )

ρ k [w(p j ) -w(p k )] + ,
where x + := max(0, x) for x ∈ R.

3) and those seeing Provider k first, who agreed to pay p k but were rejected because of Provider k's limited capacity, and who also agree to pay p j , for a total demand min 1,

w(p j ) w(p k ) ρ sp k ,
where ρ sp k is the part of the demand w(p k )ρ k that is spilled over by Provider k. The total demand ρ T j (p j , p k ) for Provider j then equals the sum of the aforementioned components:

ρ T j (p j , p k ) := w(p j )ρ j + ρ k [w(p j ) -w(p k )] + + min 1, w(pj ) w(p k ) ρ sp k
Note the dependance in both prices, although for simplicity we will sometimes just write ρ T j when there is no ambiguity. When the total demand at an RSU exceeds its capacity, some requests are rejected: we assume the RSU serves users up to its capacity level, and the rejected requests are selected randomly among all requests. This leads to an identical probability of success P j for each request submitted to Provider j, that is simply given by

P j = min 1, c ρ T j (1)
so that the served traffic at RSU j equals ρ T j P j = min(c, ρ T j ). Again, the probability P j depends on the price vector (p i , p j ). The corresponding revenue of Provider j is then

R j = p j min[c, ρ T j (p j , p k )]. (2)
The traffic ρ sp j spilled over by Provider j (and that will then enter the competitor's coverage area) also depends on both prices through the probability P j , and equals

ρ sp j = w(p j )ρ j (1 -P j ), (3)
with

P j = min   1, c w(p j )ρ j +[w(p j )-w(p k )] + ρ k + min[1, w(pj ) w(p k ) ]ρ sp k   .
Remark that for a given price configuration (p 1 , p 2 ), the success probabilities P 1 and P 2 are the solution of a fixedpoint system, since the success probability P j of Provider j depends on the spillover demand ρ sp k and thus on P k , that itself depends on ρ sp j and hence on P j . More specifically, assuming without loss of generality that p 1 ≥ p 2 , those success probabilities should satisfy

   P 1 = min 1, c w(p1)(ρ1+ρ2)-w(p1)ρ2P2 P 2 = min 1, c w(p2)(ρ1+ρ2)-w(p1)ρ1P1 . (4)
Proposition 1: For any price vector (p 1 , p 2 ), the system (4) has a unique solution (P 1 , P 2 ).

Proof: We again assume without loss of generality that p 1 ≥ p 2 . Since the right-hand sides of the equations in (4) are continuous in (P 1 , P 2 ) and fall in the interval [0, 1], Brouwer's fixed-point theorem [START_REF] Kakutani | A generalization of Brouwer's fixed point theorem[END_REF] guarantees the existence of a solution to the system.

To establish uniqueness, remark that P 2 is uniquely defined by P 1 through the second equation in (4), so

(P 1 , P 2 ) is unique if P 1 is unique. But P 1 is a solution in [0, 1] of the fixed-point equation x = g(x) with g(x) := min   1, 1 a + b -b min 1, 1 a+b+ -ax   , where a = w(p1)ρ1 c , b = w(p1)ρ2 c
, and = (w(p2)-w(p1))(ρ1+ρ2) c are all positive constants; we also assume a > 0 and b > 0 otherwise the problem is trivial. As a combination of two functions for the form x → min 1, 1 K1-K2x , g is continuous, nondecreasing, strictly increasing only on an interval [0, x] (if any) -it is in addition convex on that interval-, and constant for x ≥ x (note we can have x = 0 or x ≥ 1).

Assume g(x) = x has a solution x ∈ (0, x]. Then g is left-differentiable at x, and

g (x) = x2 ab (a + b + -ax) 2 ≤ x2 a (a + b + -ax) (5)
where we used the fact that x ≤ 1 (as a fixed point of g). Moreover, since x is in the domain where g is strictly increasing we have η := Remark also that g (x) < 1 if x < 1. We finally use the fact that g(0) > 0 to conclude that the curve y = g(x) cannot meet the diagonal y = x more than once: assume two intersection points x1 < x2 , then g (x 1 ) < 1 thus the curves cross at x1 , another intersection point x would imply g (x 2 ) > 1 (recall g is convex when strictly increasing), a contradiction. Hence the uniqueness of the fixed point and of the solution to (4). We can also establish continuity properties for the solution of (4), which will be used in the remainder of this paper.

Proposition 2: The success probability pair (P 1 , P 2 ) is continuous in the price profile (p 1 , p 2 ).

Proof: For a given price profile (p 1 , p 2 ), the solution (P 1 , P 2 ) of (4) can also be seen as a solution of the minimization problem

min (P1,P2)∈[0,1] 2 P 1 -min 1, c w(p 1 )(ρ 1 +ρ 2 ) -w(p 1 )ρ 2 P 2 2 + P 2 -min 1, c w(p 2 )(ρ 1 +ρ 2 ) -w(p 1 )ρ 1 P 1 2 ,
where the objective function is jointly continuous in (P 1 , P 2 ) and (p 1 , p 2 ). From the Theorem of the Maximum [START_REF] Berge | Topological spaces[END_REF], the mapping of prices (p 1 , p 2 ) into the corresponding set of solutions (P 1 (p 1 , p 2 ), P 2 (p 1 , p 2 )) is an upper hemicontinuous correspondence. From the uniqueness result above, that correspondence is single-valued and hence continuous. We therefore have continuity for p 1 ≥ p 2 and for p 2 ≥ p 1 (exchanging the roles of providers), hence continuity for all price profiles.

IV. REVENUE-MAXIMIZING PRICE FOR A PROVIDER

In this section we assume that provider k has already chosen his price, while provider j has to set his. We describe the revenue function of provider j for different scenarios, and provide an example when the willingness-to-pay function is linear.

In this whole section, we only consider prices p such that w(p) > 0, since a larger price would yield no revenue to the provider setting it.

We first establish a monotonicity result, that will be useful in the rest of the analysis.

Lemma 1: The total demand ρ T j of provider j is a continuous function of his price p j ; that function is in addition nonincreasing while provider j is not saturated (i.e., while ρ T j < c). Proof: Recall that

ρ T j (p j , p k ) = w(p j )ρ j + ρ k [w(p j ) -w(p k )] + + min (w(p k ), w(p j )) ρ k (1 -P k ).
The components of the first line are trivially continuous and non-increasing in p j with our assumptions on w(•).

The continuity of ρ T j (p j , p k ) follows from the continuity of P k in the price vector (p j , p k ), established in the previous section.

To establish the monotonicity result, we distinguish two cases.

• If p j ≤ p k , we show that the success probability P k is nondecreasing in p j : applying System (4) (with k = 1, j = 2) we get that P k is the solution of the fixed-point equation x = g(x), where the function g can be written as

g(x) = min    1, c w(p k )ρ k +w(p k )ρ j 1- c w(pj )(ρj +ρ k )-w(p k )ρ k x +   .
We then remark that, all else being equal, g(x) is nondecreasing in p j , so the solution P k of the fixed-point equation g(x) = x is also non-decreasing in p j .

As a result, when p k ≥ p j the component min (w(p k ), w(p j )) ρ k (1 -P k ) decreases with p j , and so does ρ T j . • If p k < p j , then we have

ρ T j (p j , p k ) = w(p j )ρ j + w(p j )ρ k (1 -P k ).
When ρ T k < c, then P k = 1 and ρ T j is non-increasing in p j . Now if ρ T k > c then from System (4) (this time with k = 2, j = 1), we have w(p k )(ρ j + ρ k ) -w(p j )ρ j P j > c and

ρ T j (p j , p k ) = w(p j )(ρ j + ρ k ) +w(p j )ρ k c w(p k )(ρ j + ρ k ) -w(p j )ρ j P j .
Assuming that provider j is not saturated, P j = 1 and thus

ρ T j = f (w(p j )) with f (x) := x(ρ j + ρ k ) -xρ k c w(p k )(ρ j + ρ k ) -xρ j .
But f is a non-decreasing function of x when x ∈ [0, w(p k )] and w(p k )(ρ j + ρ k ) -xρ j > c: differentiating we indeed get

f (x) ρ j + ρ k = 1 -ρ k c w(p k ) (w(p k )(ρ j + ρ k ) -xρ j ) 2 ≥ 1 - ρ k w(p k ) w(p k )(ρ j + ρ k ) -xρ j ≥ 1 - ρ k w(p k ) w(p k )(ρ j + ρ k ) -w(p k )ρ j ≥ 0,
where we used w(p k )(ρ j + ρ k ) -xρ j > c in the second line, and x ≤ w(p k ) in the last one. The non-increasingness of ρ T j = f (w(p j )) in p j then comes from that of w(•).

A. Capacity saturation price

For further analysis, we define the capacity saturation price of a provider, that depends on the price of his competitor.

Definition 1: The capacity saturation price of provider j is

p c j (p k ) := inf{p ∈ [0, p max ] : ρ T j (p, p k ) < c}.
Since ρ T j (p max , p k ) = 0, for all p k we have p c j (p k ) < p max . Additionally, Lemma 1 implies that if p c j > 0, then ρ T j (p c j , p k ) = c and p j ≤ p c j ⇒ ρ T j ≥ c. We now provide analytical expressions for that price, in the case when ρ

T j (p, p k ) ≥ c. In that case ρ T j (p c j ) = c, hence p c j satisfies    w(p j )ρ j + ρ k [w(p j ) -w(p k )] + + min 1, w(pj ) w(p k ) ρ sp k = c, ρ sp k = w(p k )ρ k [w(p k )-w(pj )] + ρj +w(p k )ρ k -c [w(p k )-w(pj )] + ρj +w(p k )ρ k + . (6)
Let us define a generalized inverse of w, as

W (q) := inf{p ∈ [0, p max ] : w(p) < q}. (7)
For q ≤ 1, W (q) is the maximum price that can be accepted by a proportion q of users.

Then the capacity saturation price can be computed as follows. (The proof is omitted due to space constraints.)

• If w(p k ) ≤ min[ c ρj , c ρ k ], then p c j = W c+w(p k )ρ k ρj +ρ k . • If c ρ k < w(p k ) ≤ 2c ρj +ρ k , then p c j = W 2c ρj +ρ k . • If c ρj < w(p k ) ≤ 2c ρj +ρ k , then p c j = W c ρj . • If w(p k ) > 2c ρj +ρ k , then p c j = W (x), with x the unique solution in [0, w(p k )] of -x 2 ρ j + x w(p k )(ρ j + ρ k ) -c ρ k -ρ j ρ j + ρ k -cw(p k ) = 0.

B. Piece-wise expression of the revenue function

The revenue function of each provider j is continuous in his price (from the continuity of ρ T j and of P j ), and can be expressed analytically on different segments.

1) When p j ≤ p c j (p k ) (or ρ T j (p j ) ≥ c when p c j (p k ) > 0), the RSU capacity of provider j is saturated, and thus his revenue is simply

R j = p j c. ( 8
)
This is the case in Figure 3 for prices p j below approximately 2.5. Figure 5 shows that for these prices, provider j spills some flow over toward provider k. Above p c j , provider j is not saturated anymore. Then if the total demand ρ k (p c j , p k ) of the competitor is strictly below c, we have a price range with no provider being saturated. In that case we have no spillover demand, and the revenue of provider j is:

R j = p j w(p j )ρ j + [w(p j ) -w(p k )] + ρ k .
If p c j < p k , then we remark that necessarily ρ T k (p k , p k ) ≤ c, i.e., we meet the price of the opponent provider before he gets saturated. Indeed, at (p c j , p k ) provider j does not spill traffic over to k, thus ρ T k (p c j , p k ) = ρ k w(p k ) (where we also used the fact that w(p k ) ≤ w(p c j )). From the definition of p c j , provider j is not saturated at (p k , p k ),

so that ρ T k (p k , p k ) = ρ k w(p k ) = ρ T k (p c j , p k ) ≤ c. Summarizing, we then have the two following segments. 2) If p c j < p k and ρ T k (p c j , p k ) ≤ c, then for p j ∈ [p c j , p k ] R j = p j (w(p j )(ρ j + ρ k ) -w(p k )ρ k ) .
Remark that this segment is empty if p c j ≥ p k or ρ T k (p c j , p k ) ≥ c. Figure 5 illustrates that when p j is between approximately 2.5 and 4, provider j serves his own traffic and the one from the competitor who refused the price p k but agrees to pay p j .

3) If ρ T k (p c j , p k ) ≤ c, then for p j ≥ max(p c j , p k ) we have while provider k remains unsaturated:

R j = p j w(p j ). 4) Now if ρ k (p c j , p k ) > c, then provider k is saturated for p j ∈ [p c
j , p max ] (which is easy to see since j has no spillover traffic), and for p j ∈ [p c j , p k ] we have

R j = p j (w(p j )(ρ j + ρ k ) -c) .
Remark that this segment appears only when both providers can be simultaneously saturated, a case not occurring in the example we display here.

5) There may be a price of provider j larger than p k , and above which the competitor gets saturated, so that provider j may serve part of the traffic spilled over by k. In that case the revenue of provider j is:

R j = p j w(p j )ρ j + w(p j ) w(p k ) ρ sp k , (9)
where

ρ sp k = w(p k )ρ (w(p k ) -w(p j ))ρ j + w(p k )ρ k + ρ sp j -c (w(p k ) -w(p j ))ρ j + w(p k )ρ k + ρ sp j .
(10) Figure 6 shows that provider k gets saturated, and the spillover traffic is served partly by provider j as illustrated in Figure 5. Figure 2 illustrates those different zones for the special case ρ 1 = ρ 2 = 11, c = 10, and w(p) = 1 -p/10. Figure 3 shows the corresponding different segments for R j when p k = 4, and Figure 4 for various prices p k of the competitor. We observe that a revenue-maximizing price can belong to different segments, depending on the competitor's price.

V. PROVIDERS PRICING GAME

In this section we consider a non-cooperative game, where providers -the players-simultaneously choose their prices, trying to maximize their individual payoffs given by (2). Our aim is to find a Nash equilibrium (NE) of this game: a pair of prices (p 1 , p2 ), such that no player could increase his revenue by unilaterally changing his price. Further, we investigate the situation where providers would decide to cooperate, trying to maximize the sum of their individual revenues (as a monopoly would do). We analyze how much the providers lose in terms of total revenue by refusing to cooperate.

Below is a more formal definition of the Nash equilibrium in the pricing game. Definition 2: A pair of prices (p 1 , p2 ) is a Nash equilibrium for the pricing game if

R 1 (p 1 , p2 ) ≥ R 1 (p 1 , p2 ) for all p 1 ∈ (0, p max ], R 2 (p 1 , p2 ) ≥ R 2 (p 1 , p 2 ) for all p 2 ∈ (0, p max ].
Nash equilibria can be interpreted as predictions for the outcome of the competition between selfish entities, assumed rational and taking decisions simultaneously.

A. The case of large capacities

We first consider here that RSUs capacities are larger than the users flows (c ≥ ρ j + ρ k ). So, for any price pair RSUs capacities are not saturated and spillover traffic never appears.

Without loss of generality we consider that ρ 1 = αρ 2 = αρ, for α ∈ (0, 1]. (The case α = 0 is trivial and not considered

Flow

Own SP-j Refused Fig. 6. Flow served by provider k for p k = 4. "Own" denotes the part of the original flow ρ k served by provider k, "SP-j" the part from users agreeing to pay p j but unserved by j due to capacity constraints, and "Refused" is the part of users refusing to pay p j .

here.) In all this subsection, we consider a linear willingnessto-pay function, i.e., w(p) = 1 -p/p max for some p max > 0.

We consider two cases separately: 1) When p 1 ≤ p 2 , the provider revenue functions are

R 1 = p 1 (w(p 1 )ρ(1 + α) -w(p 2 )ρ), R 2 = p 2 w(p 2 )ρ.
After getting derivatives and equating them to zero, we get the following pair of prices, as a Nash equilibrium candidate: p1 = pmax(α+1/2)

2(1+α)

, p2 = pmax 2 . Note that p2 > p1 for all α ∈ (0, 1]. As can be checked (see Appendix A), this pair of prices is indeed a Nash equilibrium for all α ∈ (0, 1]. The corresponding total revenue is

R = R 1 + R 2 =
p max ρ(α 2 + 2α + 5/4) 4(1 + α) .

2) When p 1 > p 2 , the revenue functions are:

R 1 = p 1 w(p 1 )αρ, R 2 = p 2 (w(p 2 )ρ(1 + α) -w(p 1 )αρ),
giving the NE candidate p1 = pmax 2 , p2 = pmax(1+1/2α)

2(1+α)

.

Again, for all α ∈ (0, 1], p1 > p2 is verified. But this pair of prices is a Nash equilibrium only for α ∈ [s, 1], where s ≈ 0.73 (the details are provided in Appendix A). The corresponding total revenue is R = p max ρ(5/4α 2 + 2α + 1) 4(1 + α) .

Summarizing, with large capacities the pricing game has

1 equilibrium if α ∈ (0, s), 2 equilibria if α ∈ [s, 1]
, with s ≈ 0.73.

We now compare the minimum total revenue in the duopoly case with the revenue a monopolist would obtain, to evaluate the cost of competition. Following the literature on the Price of Anarchy [START_REF] Koutsoupias | Worst-case equilibria[END_REF], we use the ratio between the total revenue in the worst-case Nash equilibrium and the monopoly total revenue as the cost measure.

It is easy to check, that the second Nash equilibria highlighted before -corresponding to the case p 1 > p 2 -gives a lower total revenue if it exists. Our cost of competition, as a function of α, therefore has two segments:

4(1+α 3 ) (3+4α)(α 2 +2α+5/4) if α ∈ (0, s), 4(1+α 3 ) (3+4α)(5/4α 2 +2α+1) if α ∈ [s, 1]. (11)
Remark that if we consider only the best-case Nash equilibrium (under a Price of Stability logic), then the first expression above applies for α ∈ [0, 1]. Figure 8 shows the cost of competition of [START_REF] Lee | Dissemination and harvesting of urban data using vehicular sensing platforms[END_REF], that is maximum for α = s, i.e., when the second candidate becomes actually an equilibrium. For computation details see Appendix B.

B. Homogeneous flows and arbitrary capacities

With arbitrary capacities, the model becomes intractable analytically. We treat here the special case when user flows are homogeneous, i.e., ρ 1 = ρ 2 := ρ. In that case, we can prove necessary conditions for a price profile to be an equilibrium.

Proposition 3: If (p j , pk ) is an equilibrium, then

pj > p c j (p k ), pk > p c k (p j ).
Proof: We first prove that if at least one provider -say jcharges a price lower than or equal to his capacity saturation price, then the price profile is not an equilibrium. Assume that (p j , pk ) is an equilibrium, with pj < p c j (p k ): then provider j is saturated and gets revenue pj c. But deviating to p c j (p k ) would improve his revenue to R j = p c j (p k )c, a contradiction. Now we prove that there is no equilibrium where at least one provider charges his exact capacity saturation price. Again we assume that (p j , p c k (p j )) is an equilibrium. From the result above we necessarily have pj ≥ p c j (p c k (p j )), hence ρ sp j = 0. • We first show that pj ≥ p c k (p j ): if it were not the case (p j < p c k , omitting writing pj in the saturation price of k)

then w(p j ) ≥ w(p c k (p j )), yielding ρ T k (p j , p c k ) = w(p c k )ρ = c, ρ T j (p j , p c k ) = 2w(p j )ρ -w(p c k )ρ = 2w(p j )ρ -c ≤ c. This implies w(p j ) ≤ w(p c k ), therefore w(p j ) = w(p c k ) thus ρ T j (p j , p c k ) = c, yielding pj ≤ p c j (p c k (p j )).
Since the opposite inequality also holds we have pj = p c j (p c k (p j )), i.e., each provider charges his saturation price. We then deduce that they are equal, because they are both maximum prices such that w(p)ρ = c, which contradicts our assumption that pj < p c k (p j ). • Therefore pj ≥ p c k (p j ). Let us now consider some p > pj . Then we have:

R j (p, p c k (p j )) = pw(p j )ρ 2 - c 2w(p c k )ρ -w(p j )ρ
.

We now prove that this revenue, as a function of p, has a positive right-derivative at p = pj . Differentiating, we get

R j (p, p c k (p j )) = (pw (p)ρ + w(p)ρ) 2 - c 2w(p c k )ρ -w(p)ρ -pw(p)ρ cw (p)ρ (2w(p c k )ρ -w(p)ρ) 2 . At (p j , p c k (p j )) the flow of provider k equals c: ρ T k (p j , p c k (p j )) = w(p c k )ρ + (w(p c k )ρ -w(p j ))ρ) = c, which implies R j (p j , p c k (p j )) = w(p j )ρ + ρp j w (p j )(1 -w(p j )ρ/c). (12)
If w (p j ) = 0, R j (p j , p c k (p j )) is strictly positive. We now show it is also the case if w (p j ) < 0.

• First, if pj > p c j (p c k (p j )), then R j (p j , p c k (p j )) = w(p j )ρ + ρp j w (p j )(1 -w(p j )ρ/c) > w(p j )ρ + ρp j w (p j ).
But as an equilibrium price pj , should maximize the revenue of provider j over (p c j (p c k (p j )), p max ), and thus pj should make the derivative of R j = pw(p)ρ equal to zero, giving w(p j )ρ + ρp j w (p j ) = 0, which implies R j (p j , p c k (p j )) > 0 . • Second, if pj = p c k (p j ), then w(p j )ρ = c and the revenue function derivative in ( 12) is equal to c > 0. Thus by increasing his price provider j could increase his revenue, therefore (p j , p c k (p j )) is not an equilibrium.

Formally, in order to show that a pair of prices is an equilibrium, we have to compare the revenue they yield with the maximum revenues in all other zones (as defined in Subsection IV-B) for each provider. Proposition 3 reduces this search, to zones where both prices are strictly above capacity saturation prices. It can be easily checked that situations where providers charge equal prices cannot be equilibria. Therefore the equilibrium candidates remaining can be characterized by

     p 1 = p 2 , p 1 > p c 1 (p 2 ); p 2 > p c 2 (p 1 ), R 1 (p 1 , p 2 ) = 0; R 2 (p 1 , p 2 ) = 0,
To show that such pairs are indeed Nash equilibria, we have to compare the revenue they give with the maximum revenue in other segments.

For a linear willingness-to-pay function, the system above only leaves two candidates

(p 1 , p2 ) ∈ {(1/2p max , 3/8p max ), (3/8p max , 1/2p max )} (13)
Numerically, we found that these two pairs of prices are indeed equilibria only when ρ/c ≤ t, with t ≈ 1.23. Figure 7 shows the best-response curves when ρ/c = t.

C. The cost of ignoring competition

In our scenario, providers may not be aware of the presence of each other (especially if they are located far from each other), and thus do not play a noncooperative game on prices. In that case each provider would treat users seeing him first the same way as he treats those coming from the competitor's direction. We estimate here the cost of this ignorance in terms of revenue loss.

Then each provider j would believe his total flow to be ρ 1 + ρ 2 independently of p j , and would therefore simply select his price so as to maximize p j min(c, (ρ 1 +ρ 2 )w(p j )), leading to a situation where p 1 = p 2 = arg max p p min(c, (ρ 1 + ρ 2 )w(p)). But from that situation, one provider could decrease his price and also serve users from the competitor who found him too expensive. For a linear willingness-to-pay function, we found that in the symmetric-flow case, this price shift would yield a 12.5% revenue increase, and in the large capacity case with heterogeneous flows, an approximate 33% revenue increase for α = 1/2 (see Appendix C).

VI. CONCLUDING REMARKS

In this work, we have analyzed the price competition between roadside units operators which are providing wireless access to moving vehicles. In contrast to general work on pricing in wireless networks where wireless operators compete on price to attract "common" users, in the reference scenario competition arises due to vehicles mobility even if the roadside units do not have overlapped coverage areas. The strategic interaction between two roadside units operators is analyzed through game theoretic tools. Namely, the analysis of the best-response function having fixed the competitor's price sheds light on interesting and counterintuitive behaviors in the systems which lead one roadside unit operator to increase her price to induce higher interference onto the competitor, thus causing traffic spillover from the competitor's side. The results from the best-response analysis are then leveraged to characterize the simultaneous competitive game between the two roadside units operators in terms of equilibrium existence and optimality. A natural follow-up for the present work includes the analysis of the scenario in which the strategy space of each roadside unit operator also includes the position of the roadside unit besides setting the price.

Fig. 1 .

 1 Fig.1. Flows involved in the model: among the total potential demand ρ j seeing Provider j first, we distinguish ρ sp j (demand from users agreeing to pay p j , but not served by this provider), ρ ref j (demand from users refusing to pay p j ) and ρ us j (demand unserved by any provider).

1 a+b+ -ax ≤ 1

 11 on one hand, and x = 1 a+b-bη on the other side. Their combination yields x ≤ 1 a and finally g (x) ≤ x ≤ 1.

10 PFig. 2 .

 102 Fig. 2.Capacity saturation prices, and the different zones where the expressions of revenues vary.

Fig. 3 .p k = 8 Fig. 4 .

 384 Fig. 3. Revenue of provider j as a function of his price p j , for p k = 4, ρ 1 = ρ 2 = 11, c = 10 and w(p) = 10-p 10 , illustrating the different segments.

Fig. 5 .

 5 Fig.5. Flow served by provider j, for p k = 4. "Own" denotes the part of original flow ρ j served by provider j, "SP-k" the part from users who agreed to pay p k but were unserved by k due to capacity constraints, and "Refused" the part from users who refused to pay p k

2 p 1 =BR 1 (p 2 ) p 2 =BR 2 (p 1 )Fig. 7 .Fig. 8 .

 211222178 Fig.7. Best-response curves in general case, when pmax = 10, for the maximum value of ρ/c such that an equilibrium exist. At the equilibrium point (3.75, 5), the best-response function of Provider 2 is discontinuity: that provider is indifferent between the maximum in the segment where Provider 1 is saturated (≈ 6.2) and the maximum in the segment where provider 1 is not saturated (3.75).

Note that we ignore vehicle-to-vehicle communications in this paper.

APPENDIX A PROOF OF THE EXISTENCE OF AN EQUILIBRIUM WITH LARGE CAPACITIES

To show that a pair of prices is a Nash equilibrium, we verify that no provider can increase his individual revenue by unilaterally changing his price. For the large capacities case, the revenue curves of both providers have only two different expressions, since RSUs are never saturated. Again, we consider two cases:

• First, when p 1 ≤ p 2 we have to verify that the price profile p1 = pmax(α+1/2) 2(1+α) , p2 = pmax 2 is an equilibrium.

For Provider 1, we thus have to check whether the revenuemaximizing price p 1 = p max /2 in the zone where p 1 ≤ p 2 could be better for him, or equivalently, determine for which α we could have R 1 (p 1 , p2 ) > R 1 (p 1 , p2 ), in which case the proposed profile is not an equilibrium. Taking the mathematical expressions:

We then observe that for any α ≥ 0, we have

Hence Provider 1 cannot increase his revenue by unilaterally changing his price.

For Provider 2, we similarly have to check, whether the peak price p 2 = pmax+ p1α 2(1+α) in the zone where p 1 ≤ p 2 could be preferable. Let us first rewrite

The revenue corresponding to this price is

while with p2 Provider 2 gets a revenue

Therefore the condition

For α ∈ [0, 1] that last condition is never satisfied (the function being nonpositive), and thus Provider 2 cannot increase his revenue by unilaterally changing his price.

As a result, the proposed price profile (p 1 , p2 ) is indeed a Nash equilibrium for all α ∈ (0, 1].

• We now turn our attention to the case when p 1 > p 2 . As before, we have to verify whether

, is a Nash equilibrium. Let us study whether Provider 1 has an incentive to change his price to p 1 = αpmax+p 2 2(1+α ) (his best price in the zone where he is cheaper than the competitor). Plugging the expression of p2 we have:

which would earn Provider 1 a revenue

while under the proposed price profile he gets

The polynomial expression above has a unique root s ≈ 0.73: thus, for α < s Provider 1 could increase his revenue by changing his price, and (p 1 , p2 ) is not an equilibrium. On the other hand, for α > s, the price p1 is a best-response of Provider 1 to p2 . Finally, for Provider 2 we follow the same logic, investigating whether taking the optimal price above the price of Provider 1 could lead to a revenue increase:

Observing that

APPENDIX B COMPUTING THE COST OF COMPETITION (FOR PROVIDERS)

We first derive an expression for the optimal (for providers) revenue value, which is the maximum possible sum of their revenues (that we can reach by collaborating). We use a linear expression for the willingness-to-pay function w(p) = pmax-p pmax As previously, we consider two cases:

• First, when p 1 ≤ p 2 the total revenue is

The necessary extremum condition are :

Therefore the prices maximizing the total revenue should satisfy

.

Solving that system gives

It is easy to verify that the sufficient conditions for this pair to be a maximum are also satisfied. The maximum total revenue with p 1 ≤ p 2 is therefore

• Let us now consider the case p 1 > p 2 . The total revenue is:

The necessary extremum conditions are :

, from which we get

Again, the sufficient maximality conditions are satisfied; the total revenue is:

Comparing both cases, we find that for α ∈ (0, 1] the total revenue with p 2 ≥ p 1 is larger than (or equal to) in the other case, meaning that the price profile (15) maximizes total revenue.

Then, we divide that revenue by the minimum equilibrium revenue (i.e., we compute the Price of Anarchy 11 for the game played among providers); it is easy to remark that the equilibrium yielding the smallest total revenue is the one with p 1 > p 2 (which exists only for α ∈ [s, 1], otherwise there is only one equilibrium).

APPENDIX C THE COST OF IGNORING THE COMPETITOR'S PRESENCE

Let us consider a situation where each provider believes he is the only one serving users, so that he can set his price to the monopoly price (p max /2 for linear willingness-to-pay functions, if capacities are sufficiently large). Then in practice each provider will serve only some of the users seeing him first. But from such a situation, one provider could lower his price to serve some of the traffic that refused to pay the price of the competitor, and increase his revenue. We compute here the amount of extra revenue that a provider could get by making this price move.

When p j ≤ p k , the revenue of Provider j with symmetric flows (ρ

The optimal price when the willingness-to-pay function w(•) is linear equals 3/8p max , leading to R j = 3/8p max (5/4ρ -1/2ρ) = 9/32p max ρ, while the revenue was R j = 1/4p max ρ initially. Hence Provider j can improve his revenue by a factor 1/8 ≈ 12%.

Let us consider again the large capacity case, but heterogeneous flows. Assume ρ 1 = α ρ < ρ 2 = ρ. When Provider 1 sets his price below p 2 = p max /2, his revenue is to be compared to R 1 = 1/4ραp max when ignoring Provider 2 (i.e., when taking p 1 = p max /2). The ratio between the two revenues equals (α+1/2) 2 α(1+α) and is maximized for α = 1/2, in which case it equals 4/3.


S Gueorguieva

R Synave
email: synave@labri.fr

Trochlear surface reconstruction and evaluation based on laser scanner acquisition

The purpose of this study is to investigate the accuracy, reproducibility and validation of linear measurements on digital reconstruction of the trochlear articular surface. Surface reconstruction is produced by a laser scanner acquisition pipeline that preserves rough data. Arc and chord measurements between chosen landmarks on physical specimens are simulated by geodesic path length evaluation on the corresponding digital models. The osteological specimens are issued from an archaeological series emerging from the "Soeurs Grises"'s cemetery (medieval sample, XV-XVIII century) located in Beauvais (Oise, France).

I. 3D IMAGING IN DIGITAL ANTHROPOLOGY RESEARCH

Medical imaging application in digital anthropology is an area of growing researches [START_REF] Aiello | Laser scanning and paleoanthropology: an example from olduvai gorge,tanzania[END_REF], [START_REF] Zollikofer | Computer-assisted paleoanthropology[END_REF], [START_REF] Recheis | New methods and techniques in anthropology[END_REF], [START_REF] Spoor | Using diagnostic radiology in human evolutionary studies[END_REF], [START_REF] Subsol | 3d image processing for the study of the evolution of the shape of the human skull: presentation of the tools and preliminary results[END_REF], [START_REF] Weber | Virtual anthropology (va): a call for glastnost in paleoanthropology[END_REF], [START_REF] Zonneveld | Applications and pitfalls of ct-based 3d imaging of hominid fossils[END_REF], [START_REF] Mafart | Three-dimensional computer imaging of hominid fossils: a new step in human evolution studies[END_REF]. High resolution 3D digital representations enable accurate reconstruction and geometric morphometrics analysis of specimens. Computed tomography is presently acknowledged as one of the most important technique in digital anthropology. In particular, the reconstruction and the restoration of human fossil is widely viewed as a tool for cultural heritage preservation [START_REF] Spoor | Using diagnostic radiology in human evolutionary studies[END_REF], [START_REF] Bruner | Digital tools for the preservation of human fossil heritage: Ceprano, saccopastore, and other case studies[END_REF]. Often, virtual reconstruction stands for the visualization of complex inner anatomical structures as cortical bone thickness [START_REF] Spoor | Linear measurements of cortical bone and dental enamel by computed tomography: applications and problems[END_REF], inner ear and endocranial volumes [START_REF] Zollikofer | Computer-assisted paleoanthropology[END_REF], [START_REF] Ayache | Epidaure: A research project in medical image analysis, simulation and robotics at inria[END_REF]. Further, rapid prototyping produces 3D solid reproductions of the anatomical models [START_REF] Mankovich | The display of threedimensional anatomy with stereolithographic models[END_REF], [START_REF] Seidler | A comparative study of stereolithographically modelled skulls of petralona and broken hilll: implications for future studies of middle pleistocene hominid evolution[END_REF], [START_REF] Flaczynski | Modélisation informatique et stéréolithographique d'un crâne sec nubien daté de 2700 ans et présentant une craniosténose[END_REF], [START_REF] Zelicourt | Single-step stereolothography of complex anatomical models for optical flow measurements[END_REF], [START_REF] Allard | Use of hand-held laser scanning and 3d printing for creation of a museum exhibit[END_REF]. The acquisition, reconstruction and reproduction are wellestablished as qualitative techniques. Less is done in evaluating of their potential as "quantitative tool". A great number of researches [START_REF] Jungers | Computed tomography and biomechanical analysis of fossil long bones[END_REF], [START_REF] Tate | High-resolution computed tomography for the comparative study of fossil and extant bone[END_REF], [START_REF] Sumner | Computed tomography and automated image analysis of prehistoric femora[END_REF], [START_REF] Ruff | Use of computed tomography in skeletal structure research[END_REF], [START_REF] Richtsmeier | Precision, repeatability, and validation of the localization of cranial landmarks using computed tomography scans[END_REF], [START_REF] Hildebolt | Validation study of skull three-dimensional computerized tomography measurements[END_REF], [START_REF] Spoor | Linear measurements of cortical bone and dental enamel by computed tomography: applications and problems[END_REF] illustrates the complexity of the measurement definition with respect to its precision, repeatability and validation. As for example, Spoor's work [START_REF] Spoor | Linear measurements of cortical bone and dental enamel by computed tomography: applications and problems[END_REF] assesses linear measurements of cortical bone thickness by computed tomography. The cortical bone CT scan measurements are compared with the measurements on the corresponding cortical bone physical sections in order to validate the range values. The reported maximum error range is of ±0.1mm and fit to the "best possible accuracy" depending on the CT characteristics and reconstruction settings in use (slice thickness of 1.5mm, spatial resolution is within the interval [0.5mm, 0.6mm], and pixel size within [0.16mm, 0.22mm] (see also [START_REF] Zonneveld | The technique of direct multiplanar high resolution ct of the temporal bone[END_REF]). Moreover, CT based measurements are performed in the same slice image. Localization of the cranial landmarks using CT scans is studied by [START_REF] Richtsmeier | Precision, repeatability, and validation of the localization of cranial landmarks using computed tomography scans[END_REF]. A "gold standard" set of data, recording the relative positions of 35 landmarks on a collection of ten skulls, is used to evaluate the precision of the landmark locations. Measurements are carried out inside a single CT scan as well as inside two CT scan images. While the average difference in landmark locations within an CT scan is less than 0.5mm for the entire gold standard set, some distance measures within two CT scans have range error greater than 2.5mm. An error distortion along the coordinate axes is also emphasized. Disparity between errors along the axes is an important drawback since it could invalidate the identification of outlines defined on the basis of particular space orientation. Further, using "fuzzy landmarks" [START_REF] Valeri | Capturing data from three dimensiona surfaces using fuzzy landmarks[END_REF], the average error could be reduced to 1.15mm. More recently, a comparison between mandibular landmarks position computations from CT scans and by a 3D space digitizer [START_REF] Williams | Comparaison of mandibular landmarks from computed tomography and 3d digitizer[END_REF] establishes a standard deviation in landmark coordinates of 1.136mm with respect to linear distance range within the limits of [0.001mm, 3.889mm]. The reported and discussed above distance measurements are linear distances. Neither curvilinear arc distance measurements nor surface area measurements are experimented as long as we know. Several authors conclude that in place of slice images, 3D reconstructions could result in error reduction. No matter of the support in use, physic specimen, numerical representation or synthetic reproduction, a geometric morphometrics based analysis requires precise localization of landmarks and measurements of not only of linear distances, but also of surface area and volume data features and regions of interest. In the present work, rough data acquisition is based on the the laser scanner acquisition pipeline. Our goal is to demonstrate that laser scan based 3D reconstruction have good capacity for quantitative assessment of linear and curvilinear distances. Moreover, the mesh models that are a direct output of the acquisition pipeline can support iD, 0 ≤ i ≤ 3, axis orientation independent linear and curvilinear measures within comparable precision and repeatability. All measures are validated in comparison with anthropometric measures taken with caliper and millimeter ruler band on the specimens.

II. MATERIAL AND METHODS

LANDMARK DATA FROM ARCHAEOLOGICAL SERIES The studied tali are issued from an archaeological series emerging from the "Soeurs Grises"'s cemetery (medieval sample, XV-XVIII century) located in Beauvais (Oise, France). An almost exhaustive excavation was realized in 1992 and allowed the bring up to day of 350 inhumations and of about fifty ossuaries and reduced body individuals. Up to day, the specimen sample is composed of 350 adult specimens. In the current study our interest is centered on a particular concavo-convex (trochlear) articular facet on the superior surface of the talus's body. In almost all specimens (10 pairs of left and right tali) the trochlea is relatively well preserved. It is strongly convex in the anterior direction, moderately concave in the mediolateral plane and wider distally than proximally. Anatomically, the trochlea is very important articular facet since it is involved in the ankle joint and consequently in the foot position during locomotion [START_REF] Wood | A technique for establishing the identity of "isolated" fossil hominin limb bones[END_REF]. The definitions and the locations of the landmark points follow [START_REF] Saller | Lehrbuch der Anthropologie[END_REF], [START_REF] Scheuer | The development of juvenile osteology[END_REF], [START_REF] De | Utilisation d'outils de mesures en 3d dans le cadre d'une etude comparative morphofonctionnelle de tali d'hominoides actuel et du hominide fossile[END_REF]. Three pairs of landmark points for the trochlear shape are defined: (P apmr , P aplr ), (P ipmr , P iplr ) and (P ppmr , P pplr ). Each pair corresponds to one position: the anterior, the intermediate and the posterior with respect to the frontal plan. A point of the pair is either on the medial or on the lateral process of the trochlea with respect to the sagittal plane, denoted as medial or lateral ridge. Landmark point location for the specimen BCT 92 -S319 -G are given in Fig. 1.

The anthropometric arc measures are with Mitutoyo digital caliper resolution of 0.01mm. The chord measures are taken with a DKSH millimeter ruler band with expected accuracy of 1mm. All the measurements quoted with "M" have been defined by Martin and Saller (1957) [START_REF] Saller | Lehrbuch der Anthropologie[END_REF] and the others are new measurements introduced for the present study. Let us define the anterior trochlea ridge as the curve starting in P apmr passing through the anterior extremity of the medial lengthwise curvature of the trochlea and finishing in P aplr . By analogy, the posterior trochlea ridge is defined between P ppmr and P pplr . The following landmark point distances are -M 5 is the length of the chord P ipmr P iplr ; -M 5(1) is the length of the chord P ppmr P pplr ; -M 5(2) is the length of the chord P apmr P aplr ; -Cm is the chord length P ppmr P apmr ; -Cl is the chord length P pplr P aplr ; -Am is the length of the arc joining P ppmr and P apmr along the medial ridge; -Al is the length of the arc joining P pplr and P aplr along the lateral ridge. The illustration is given in Fig. 1 and Fig. 2.

The error standard deviations (STD) of manual linear measurements on the series of 20 specimens are shown in Fig. 3. As it is seen only 10.83% of the total number of chord measures have STD values under 1mm, 40.83% are under 2mm, and 74.17% are under 3mm.

In the pursuit of the current knowledge, measures on two individuals, BCT 92 -S319 -G and BCT 92 -S380 -G, will be expressly stated but issues for the rest of the studies talus collection are within the same limits. Range values for manual measures on the specimens BCT 92 -S319 -G and BCT 92 -S380 -G are given in Fig. 5 and Fig. 6.

The STD of manual arc measurements on the series of 20 specimens are shown in Fig. 4. It could be seen that the curvilinear manual measures have bad performance: 5% of the total number of arc measures have STD values under 1mm, 22.5% are under 2mm, and 55% are under 3mm.

LANDMARK DATA FROM IMAGE PROCESSING

The experimental platform in use consists of a 3D digitizer VIVID 300/VI-300, a software package A2RI and a printing system EDEN250. The acquisition is done by a non-contact 3D digitizer VIVID 300/VI-300 with a laser wavelength λ = 690nm and maximal optical power P M ax = 7mW , object distance range in the limits of [0.55m, 1.2m], field of view in the limits of [185mm, 395mm] and output data points 400 × 400. In order to optimize the quality of the laser scans an acquisition protocol for the specimen acquisition is established [START_REF] Synave | Reconstruction de solides à partir d'acquisitions surfaciques[END_REF] in such a way that laser power is fixed to 15% of P M ax, object distance is 550mm, and the elongated triangles with aspect ratio superior to 4 of the scan borders are eliminated. During the acquisition, each specimen is exposed to the laser scan on a rotating stage set. The rotating stage is turned over the Z axis by a step of 45deg and thus producing 8 scans exported as triangular mesh models. The software toolkit that supports the pipeline from the acquisition to the 3D printing is the A2RI1 [START_REF] Synave | Toolkit for registration and evaluation for 3d laser scanner acquisition[END_REF]. For each specimen, the input of A2RI is the set of 8 triangular mesh models, one for each laser scan. The output is a 3D solid model of the specimen that could be exported in ST L format for the printing by the 3D printing system. It should be emphasized that all treatments during the A2RI processing do not approximate rough data. Thus all output mesh vertices correspond to acquired points on the specimen. In phantom studies [START_REF] Synave | Reconstruction de solides à partir d'acquisitions surfaciques[END_REF] it has been shown that the precision of A2RI is within the limits of [0.72mm, 1.2mm] for linear distances and of [0.68mm, 1.46mm] for curvilinear geodesic distances. Finally, the 3D printing system Eden250 used for the specimen reproduction has a "typical" accuracy in the limits of [0.1mm, 0.2mm]2 . The synthetic reproduction of BCT 92 -S319 -G processed by our pipeline is given in Fig. 1 (right). Chord measures on the numerical representations of the specimen are calculated as Euclidean distances between the vertices corresponding to the landmark points (cf Fig. 2 (left)). Arc measures fit to the length of the geodesic paths that join these vertices and that follow the surface ridges (cf Fig. 2 (right)). Range values for numerical measures on the 3D reconstructions of the specimens BCT 92 -S319 -G and BCT 92 -S380 -G are given in Fig. 7 and Fig. 8.

The STD of the measures on the 3D reconstructions of BCT 92-S319-G and BCT 92-S380-G are given in Fig. 9 and Fig. 10. For chord measures on the BCT 92 -S319 -G reconstruction, 56, 67% are under 1mm, 86, 67% are under 2mm and 96, 67% are under 3mm. For the BCT 92 -S380 -G reconstruction, 56, 67% are under 1mm, 83.33% are under 2mm and 100% are under 3mm. With respect to arc measures on the BCT 92 -S319 -G representation, 30% are under 1mm, 70% are under 2mm and 90% are under 3mm. For the The anterior width M 5(2) of BCT 92 -S380 -G joining the medial and the lateral ridges respectively in P apmr and P aplr presents a damaged neighborhood around P apmr that mislead the laser scanner acquisition and reconstruction in this region. The range intervals are comparable to the ones for the landmark point coordinate calculation given by [START_REF] Richtsmeier | Precision, repeatability, and validation of the localization of cranial landmarks using computed tomography scans[END_REF] and [START_REF] Valeri | Capturing data from three dimensiona surfaces using fuzzy landmarks[END_REF] and obtained from CT slide images.

With respect to the reproducibility, numerical measures are more stable as long as in all experiments the population of numerical measures with STD under given threshold is greater than the one for manual measures. For example, 80% of numerical measures have STD under 2mm while for the manual measures only 36, 6% achieve this precision.

The validity of the numerical measures is tested in comparison with the manual measures. Table I recapitulates results for both BCT 92-S319-G and BCT 92-S380-G specimens. A visual feedback of the 3D models and measures used in the table is given in Fig. 12. The absolute differences between numerical and manual measure values are given in the 4rd and the 7th columns. It can be seen that all absolute differences for linear distances are under 2mm except for M 5(1) on BCT 92-S380-G. In Fig. 13 from the medial and posterior view of BCT 92 -S380 -G it Arc measures for BCT 92 -S380 -G are less precise. Indeed, as it seen from the medial and lateral profiles in Fig. 11, ridge extremities are damaged that make difficult the landmark point location.

IV. DISCUSSION AND CONCLUSION

From an anatomical point of view one can conclude that the posterior width M 5( 1) is stable all over the processing and could be retained as representative trochlear surface shape feature. Further, measures on the medial trochlea ridges are less exact. The chord measures are always better than the arc ones. This implies that the representative points location is quite precise but the distances measured along trochlea ridges necessitate more accurate computation. That is important because arc distances and curvature along medial and lateral trochlea ridges are related to functional anatomy of the human foot.

For the time being, our interest is directed to shape surface estimation in regions with constraint curvature variation corresponding to contact area of the talus articular surface.

Our believe is that 3D laser scanner based image acquisition and reconstruction are well-fitted to a precise quantitative estimation of surface shape features as linear and curvilinear distances, curvature and surface area appreciation.

Fig. 1 .

 1 Fig. 1. Dorsal view of the trochlear surface of BCT92-S319-G and landmark point locations : (left) original specimen -(right) synthetic reproduction.

Fig. 2 .

 2 Fig. 2. Dorsal view of the trochlear surface of BCT92-S319-G and landmark point locations : (left) laser scanner acquisition -(right) solid mesh model.

Fig. 3 .

 3 Fig. 3. STD of chord manual measurements on the series of 20 specimens.

Fig. 4 .

 4 Fig. 4. STD of arc manual measurements on the series of 20 specimens.

Fig. 5 .

 5 Fig. 5. Chord and arc manual measures of the specimen BCT92-S319-G.

Fig. 6 .

 6 Fig. 6. Chord and arc manual measures of the specimen BCT92-S380-G.

Fig. 7 .

 7 Fig. 7. Chord and arc numerical measures on 3D reconstruction of the specimen BCT92-S319-G.

Fig. 8 .

 8 Fig. 8. Chord and arc numerical measures on 3D reconstruction of the specimen BCT92-S380-G.

Fig. 9 .

 9 Fig. 9. STD of chord and arc measures on the 3D reconstruction of BCT92-S319-G.

Fig. 10 .

 10 Fig. 10. STD of chord and arc measures on the 3D reconstruction of BCT92-S380-G.

Fig. 11 .

 11 Fig. 11. BCT92-S319-G: (left) Medial view -(right) Lateral view

Fig. 13 .

 13 Fig. 13. BCT92-S380-G: (left) Medial view -(right) Posterior view

A2RI could be download at http : //liba2ri.f ree.f r

According to the technical specifications and depending on the geometry, part orientation and print size.

ACKNOWLEDGMENT

The authors wish to thank Mr. P. Courtaud (UMR PACEA, LAPP) for providing the collection from Pessac osteological repository, and Mr. J.-M. Femolan in charge of the complete collection from Archaeological Department of Beauvais city concil.


R † And Couture-Veschambre

Ch ‡ †umr5800

Teaching geometric modeling algorithms and data structures through laser scanner acquisition pipeline

Keywords: geometric modeling, digital paleoanthropology, cultural heritage preservation, laser scanner acquisition, image registration, image reconstruction, 3D printing

Experiences from geometric modeling course based on a specific teaching medium, namely trochlear surface reconstruction from laser scans, its evaluation in terms of shape feature measurements and finally its instantiation through 3D printing, are presented. Laser scanner acquisition, reconstruction and 3D printing lend well to teaching general concepts in geometric modeling for several reasons. First, starting and ending with real physical 3D objects (the talus and its 3D print) provide in addition to the classical visual feedback a material feedback for correctness of treatments all over the pipeline. Second, the notion of error during each step of the pipeline is illustrated in a very intuitive way through length measurements, manual ones with calipers on the tali, and numerical ones with arc and chord lengths on the numerical reconstructions. Third, students are involved with challenging scientific problems and produce semester-long projects included in larger scaled project of cultural heritage preservation. Our believe is that this approach gives a deeper understanding of both theoretical and application issues in geometric modeling.

Introduction

Geometric modeling evolves into a great variety of visual computing applications as 3D shape matching and recognition [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF][START_REF] Mori | Efficient shape matching using shape contexts[END_REF], medical image analysis [START_REF] Ayache | Epidaure: A research project in medical image analysis, simulation and robotics at inria[END_REF] and cultural heritage preservation [START_REF] Levoy | The digital michelangelo project[END_REF][START_REF] Bernardini | The 3d model acquisition pipeline[END_REF][START_REF] Bernardini | Building a digital model of michelangelo's florentine pietà[END_REF][START_REF] Rusinkiewicz | Real-time 3d model acquisition[END_REF][START_REF] Coqueugniot | Range image processing for paleoanthropology heritage preservation[END_REF]. Geometric models capture spatial aspects of the objects of interest for an application [START_REF] Requicha | Representation for rigid solids: Theory, methods and systems[END_REF][START_REF] Requicha | Solid modelling and beyond[END_REF][START_REF] Requicha | Geometric Modeling: A First Course[END_REF] and necessitate rigor mathematical foundations [START_REF] Mortenson | Geometric Modeling[END_REF][START_REF] Faux | Computational Geometry for Design and Manufacture[END_REF][START_REF] Farin | Gurves and Surfaces for Computer Aided Design[END_REF][START_REF] Fomenko | Algorithms and Computer Methods for Three-Manifolds[END_REF][START_REF] Farin | Handbook of Computer Aided Design[END_REF][START_REF] Agoston | Computer graphics and geometric modeling: mathematics[END_REF][START_REF] Agoston | Computer graphics and geometric modeling: implementation and algorithms[END_REF].

Unfortunately, the mathematical prerequisites are often unappealing to college students and the visual representation of 3D objects is usually considered as a final goal. Often, facet (planar face) models correspond to "soup" of polygons with degenerate polygons, holes and self-intersecting faces. Such models could be acceptable for rough graphical representation but when looking for surface feature estimation as the surface normal for example, these models lead to erroneous calculation. A possible way to handle such irregularities is to locate artefacts through geometric intersection tests [START_REF] Coxeter | Geometry revisited[END_REF][START_REF] O'rourke | Computational geometry in C (2nd edition)[END_REF], to refine the underlying mesh in singular points and then to produce a manifold surface subdivision. These treatment is a commonly employed operation when integrating different scan views in a laser scanner acquisition based surface reconstruction. The advantage in making use of the laser scanner acquisition pipeline is that theoretical notions as rsets and 3D manif olds, topological invariants and Euler characteristics, have an imminent impact in the results. Another example of geometric model prevailing property is the validity to argue if the model corresponds to (at least) one object. In the terms of the three level hierarchical view of modeling [START_REF] Requicha | Geometric Modeling: A First Course[END_REF][START_REF] Hoffmann | Geometric and Solid Modeling: An Introduction[END_REF], "physical object → mathematical models of objects → representations of objects", a valid geometric model could correspond to a mathematical model but not to a physical object. A new intuition on validity is brought to an active state through the 3D printing. The pipeline supports a three level graph view, "physical object ↔ mathematical models of objects ↔ representations of objects". Once the object boundary surface is reconstructed, a STL file [START_REF] Szilvasi-Nagy | Analysis of stl files[END_REF] could be produced and a synthetic physical object reproduction, a 3D print, could be output. The distinction between valid and "printable" geometric model allows to close up the cycle "physical object ↔ mathematical models of objects ↔ representations of objects ↔ physical object". Each level of these modeling view is clearly identified through the pipeline and practitioner students can themselves evaluate the quality of all intermediate results and the final product. The control is done through manual measurements on the physical supports (bone specimen or its 3D print) and through numerical measurements on the numerical object reconstruction. The purpose of the present research is to show our experience in how to teach the algorithms and data structures underlying the pipeline. Our objective is to emphasize the practice interest of geometric modeling theoretical foundations for both computation implementation and anthropology investigations. Our believe is that this approach will engage students in learning geometric modeling through the development of solutions to scientific relevant problems. In the following, the main steps are explicated with the corresponding accompanying references. This program is proposed as a MS degree course of geometric modeling in the Department of Computer Science at the University of Bordeaux 1. Examples of the student project final results are also provided.

Laser scanner acquisition, registration and integration

There are several examples of 3D model acquisition pipelines [START_REF] Curless | A volumetric method for building complex models from range images[END_REF][START_REF] Levoy | The digital michelangelo project[END_REF][START_REF] Bernardini | The 3d model acquisition pipeline[END_REF][START_REF] Rusinkiewicz | Real-time 3d model acquisition[END_REF][START_REF] Jaeggli | Online 3d acquisition and model integration[END_REF][START_REF] Pauly | Example-based 3d scan completion[END_REF][START_REF] Synave | Toolkit for registration and evaluation for 3d laser scanner acquisition[END_REF]. In our case, the acquisition is done by a non-contact 3D digitizer VIVID 300/VI-300 with a laser wavelength λ = 690nm and maximal optical power P M ax = 7mW , object length range is in the limits of [0.55m, 1.2m], field of view is in the limits of [185mm, 395mm] and output data points 400 × 400. An illustration of talus acquisition with VIVID 300/VI-300 and Polygon Editing Tool is given in Fig. 1. Two views of the talus in Fig. 1(left) are registered in Fig. 1(center). It should be noted that the rough acquisition for a view corresponds to the mesh model shown in Fig. 2. The human operator is supposed to clean up those parts that are out of the region of interest and to produce the model given in Fig. 1(left). The final step is the integration of all scans and bringing forward for consideration the reconstructed surface as shown in Fig. 1(right). Acting with the system for a short time is sufficient to identify crucial points as:

-The optimal values for the laser power (% of P M ax ?) and the object distance.

A few protocols including technical data sheets for the laser scanner measurements are known [START_REF] Synave | Reconstruction de solides à partir d'acquisitions surfaciques[END_REF]. -The initial alignment and the partial overlap between sequential scans.

The widely used in range image registration Iterative Closest Point(ICP) algorithm and its derivatives [START_REF] Besl | A method for registration of 3-d shapes[END_REF][START_REF] Trucco | Robust motion and correspondence of noisy 3-d point sets with missing data[END_REF][START_REF] Rusinkiewicz | Efficient variants of the icp algorithm[END_REF][START_REF] Tsakiri | Laser scanning issues for the geometrical recording of complex statue[END_REF][START_REF] Chetverikov | Robust euclidean alignment of 3d point sets: the trimmed iterative closest point algorithm[END_REF][START_REF] Synave | Automated trimmed iterative closest point algorithm[END_REF] strongly depend on both characteristics. -The extent and the constraints to interpolate and/or approximate cloud points during the integration. Using Polygon Editing Tool for example, one can easily observe surface artefacts as undesired smoothing and lost of details in the reconstructed surface. 3 Object boundary reconstruction and evaluation

Programming projet

At this period, students start with a programming project dealing with surface reconstruction from 3D range images [START_REF] Bernardini | The ball-pivoting algorithm for surface reconstruction[END_REF][START_REF] Allen | The space of human body shapes: reconstruction and parameterization from range scans[END_REF][START_REF] Nehab | Efficiently combining positions and normals for precise 3d geometry[END_REF] obtained through the first stage. The range images could represent individual or integrated views of the studied objects and the goal is to construct valid surface mesh models. The following requirements are implied:

-The data structures in used [START_REF] Samet | Foundations of multidimensional and metric data structures[END_REF] should support:

1. Efficient topological queries as variable sized vertex neighbourhood construction referring to the quad-edge [START_REF] Guibas | Primitives for the manipulation of general subdivisions and the computation of voronoi diagrams[END_REF], the winged-edge [START_REF] Baumgart | A polyhedron representation for computer vision[END_REF][START_REF] Braid | The sythesis of solids bounded by many faces[END_REF] or the half-edge [START_REF] Mantyla | Geometric and Solid Modeling: An introduction[END_REF] data structures.

Evaluated boundary representation for rapid visualisation.

-The models should fulfil mesh element quality criteria [START_REF] Shewchuk | What is good linear element? interpolation, conditioning, and quality measures[END_REF][START_REF] Pébay | Analysis of triangle quality measures[END_REF][START_REF] Cabello | Toward quality surface meshing[END_REF][START_REF] Knupp | Remarks on mesh quality[END_REF] and geometric cohererence of data [START_REF] Botsch | A robust procedure to eliminate degenerate faces from triangle meshes[END_REF][START_REF] Bischoff | Towards robust broadcasting of geometry data[END_REF].

In order to broach the geometric computations students are urged to experiment with examples of robustness problems [START_REF] Yap | Towards exact geometric computation[END_REF][START_REF] Yap | Towards robust geometric computation[END_REF][START_REF] Yap | Special issue on robust geometric algorithms and their implementations[END_REF][START_REF] Mehlhorn | Robust Geometric Computation[END_REF][START_REF] Kettner | Classroom examples of robustness problems in geometric computations[END_REF] and in particular the Core [START_REF] Karamcheti | A core library for robust numeric and geometric computation[END_REF][START_REF] Li | Core library tutoroal[END_REF] and LEDA [START_REF] Mehlhorn | LEDA[END_REF] libraries. Along with geometric intersection tests essential in quality mesh improvement, a set of mesh query based operations are also required:

-Calculation of Euler characteristics of the underlying object boundary surface [START_REF] Alexandroff | Elementary Concepts of Topology[END_REF][START_REF] Henle | A Combinatorial Introduction to Topoly[END_REF][START_REF] Kinsey | Topology of surfaces[END_REF].

-Geodesic path construction between a pair of source and destination vertices [START_REF] Mitchell | The discrete geodesic problem[END_REF][START_REF] Kimmel | Computing geodesic paths on manifolds[END_REF][START_REF] Chen | Shortest paths on a polyhedron; part i: computing shortest paths[END_REF][START_REF] Kanai | Approximate shortest path on a polyhedral surface and its applications[END_REF][START_REF] Novotni | Computing geodesic paths on triangular meshes[END_REF][START_REF] Martinez | Computing geodesics on triangular meshes[END_REF][START_REF] Surazhsky | Fast exact and approximate geodesics on meshes[END_REF][START_REF] Polthier | Straightest geodesics on polyhedral surfaces[END_REF][START_REF] Lee | Straightest paths on meshes by cutting planes[END_REF][START_REF] Synave | Constraint shortest path computation on polyhedral surfaces[END_REF].

Finally, as an optional functionality, the discrete curvature estimation [START_REF] Polthier | Polyhedral surfaces of constant mean curvature[END_REF][START_REF] Garimella | Curvature estimation for unstructured triangulations of surfaces[END_REF][START_REF] Surazhsky | A comparison of gaussian and mean curvatures estimation methods on triangular meshes[END_REF][START_REF] Garland | Notes on differential geometry[END_REF] at mesh model vertices is required. This feature is helpful to indicate the characteristic vertex corresponding to chosen landmark as long as it is often situated either on spherical or hyperbolic point. In case when this functionality is not implemented the choice is made by visual appreciation. These programming projects are developed as a collaborative work within two or three person team. Weekly, each team present the advancement of the project. The implementation is done in C + + using either Qt4 or GT K graphical user interface [START_REF] Blanchette | C++ GUI programming with Qt4[END_REF][START_REF] Krause | Foundations of GTK+ Development[END_REF].

Validation

Once the student software becomes operating the validation of the boundary surface reconstruction starts with geometric features evaluation. First, manually on the studied specimen and then by simulating the same measurements on the numerical representation.

For our experience we study osteological specimens issued from an archaeological series emerging from the "Soeurs Grises"'s cemetery (medieval sample, XV-XVIII • century) located in Beauvais (Oise, France). In particular we are interested in the (trochlear) articular facet on the superior surface of the talus's body involved in the ankle joint and consequently in the foot position during locomotion [START_REF] Wood | A technique for establishing the identity of "isolated" fossil hominin limb bones[END_REF]. An illustration is given in Fig. 3. Three pairs of landmark points for the trochlear shape are defined following [START_REF] Saller | Lehrbuch der Anthropologie[END_REF][START_REF] Scheuer | The development of juvenile osteology[END_REF][START_REF] De La Villetanet | Utilisation d'outils de mesures en 3d dans le cadre d'une etude comparative morphofonctionnelle de tali d'hominoides actuel et du hominide fossile[END_REF]: (P apmr , P aplr ), (P ipmr , P iplr ) and (P ppmr , P pplr ). Each pair corresponds to one position: the anterior, the intermediate and the posterior with respect to the frontal plan. A point of the pair is either on the medial or on the lateral process of the trochlea with respect to the sagittal plane, denoted as medial or lateral ridge as illustrated in Fig. 4. Anthropometric chord and arc length mesurements are performed manually with caliper and millimetre ruler band on the specimens.

For the validation of the trochlear surface reconstruction each manual measurement is replicated on the mesh model. Chord measurements are similar to euclidean distances between the pair of characteristic vertices, each one corresponding to an anatomical landmark point. Arc measurements agree with geodesic distances between the characteristic vertices. Statistical analysis [START_REF] Slice | Modern morphometrics in physical anthropology[END_REF] is performed over all samples in order to evaluate the range interval for the standard deviations for the measured entities. All results are presented and appreciated by an expert in the application field namely an anthropologist. According to our experience these discussions are particularly fruitful as long as the visual and the numerical feedback of the trochlear surface reconstruction are commented by a potential final user of such kind of software package. In this way the developer is also confrontated with the technical trade specifications.

Offsetting object boundary and 3D printing

Creating of 3D print models becomes very popular as research and training tool [START_REF] Allard | Use of hand-held laser scanning and 3d printing for creation of a museum exhibit[END_REF][START_REF] Gill-Robinson | Using image analysis software to create a physical skull model for the facial reconstruction of a wrapped akhmimic mummy[END_REF][START_REF] Seulin | 3d tools for ancient wooden stamps analysis and conservation[END_REF][START_REF] Carter | The role of printing in teaching and education in human skeletal anatomy[END_REF]. The 3D printer put to student's disposal is an Objet Eden250 3D printing system. The input of the 3D printing system is the STL format [START_REF] Szilvasi-Nagy | Analysis of stl files[END_REF] that does not guaranty the correctness of the underlying surface namely to correspond to an oriented 2D manifold without boundaries. The result of the reconstruction on the previous stage is a triangulated surface but with possible existing boundaries as for example when the object boundary is partially reconstructed and does not enclose a bounded volume (shell surface). For such cases initial surface should be offset [START_REF] Rossignac | Offsetting operations in solid modelling[END_REF][START_REF] Lee | Offsetting operations on non-manifold boundary representation models with simple geometry[END_REF][START_REF] Maekawa | An overview of offset curves and surfaces[END_REF][START_REF] Chen | A point-based offsetting method of polygonal meshes[END_REF][88][START_REF] Lee | Offsetting operations on non-manifold topological models[END_REF]. The offset direction is chosen opposite to the normal direction in order to preserve the chord and arc metric of the initial shell surface. The offsetting and the STL output are added as new functionalities to the programming project. Finaly, each student team produces and evaluates the 3D print of the studied talus specimen. Evaluating the precision of the complete pipeline from the acquisition through the reconstruction to the 3D printing is an innovative approach to validate the whole treatment channel and the supported data structure and algorithm implementations. A few similar works are known related to vision system based on computed tomography and commercial software [START_REF] Nizam | Dimesional accuracy of the skull models produced by rapid prototyping technology using stereolithography apparatus[END_REF]. According to the technical specifications the typical tolerance of Eden250 is of 0.1mm. The student experimental results show that the variance of manual measurements on the 3D print is similar to the one of manual measurements on the original specimen but that the absolute values of chord and arc lengths increase for the 3D print. It is observed that the interchangeability of both supports is limited to a resolution of 1.5mm.

Discussion and experimental results

The elaborated training program is successfully applied as a graduated course in the master of Computer Science at the University of Bordeaux 1. It should be noticed that mathematical prerequisites and namely notions from algebraic topology of surfaces and differential geometry are difficult to apprehend as long as in the previous Computer Science degree courses students do not broach these subjects. Often, problems are identified but the mathematical apparatus to resolve them takes too long time to become familiar with enough in order to be used. Dealing with the laser acquisition pipeline as a whole permits to remain flexible with respect to the pedagogical objectives that can be adapted according to individual skills of the members in a programming project team. For the basic pipeline processing steps: registration, reconstruction, offsetting, geodesic path computation and curvature estimation, rough algorithms are proposed as initial guess to be improved as students go along. See for example, in Fig. 5(b) the illustration of the shortest path calculation between the source (P pabl ) and the destination (P ppbl ) characterstic vertices. The rough algorithm is the Dijkstra algorithm. Very quickly a shortest path computation could be implemented and visualized but the limits of a naive implementation are obvious even for mesh models with a total number of vertices under the range of 10 4 . Further, the problem of non-uniqueness of the shortest path and the lack of a measurement direction as in physical measurements with caliper or millimetre ruler band could be seen in Fig. 5(a). Finally, the choice of the path extremity in a specified mesh vertex position could be tedious and erroneous task so that the possibility of interaction with vertices in a variable sized neighborhood of source (destination) vertex is useful as shown in Fig. 5(c). The offsetting of a mesh model from Fig. 6(a) is illustrated in Fig. 6(b) and Fig. 6(c). The rough algorithm consists in a vertex based translational extrusion. As long as for 3D printing of shell surface we need a narrowly banded extruded surface, the offset vector is of small magnitude and in almost all tested models self-intersections are avoided. According to our experience discrete curvature evaluation stands up as one of the most apprehended problem. Starting with [START_REF] Garimella | Curvature estimation for unstructured triangulations of surfaces[END_REF] and the GSL GNU Scientific library for matrix operations, an initial evaluation could be reported as shown in Fig. 7 where flat like regions are displayed in blue color and the curved ones in different graduations of the green color. There are two check on points for the developed software: First, evaluate the reconstructed object boundary in terms of chord and arc measurements between characteristic vertices on the underlying mesh model, and second, estimate the same geometric features on the corresponding 3D print. The computed values are compared with the manual measurements on the original specimen. The goal is to establish the degree of interchangeability along the different object reproductions, the numerical representations and/or the synthetic ones. The results are summed up for the specimen the BCT92-S319-G in Fig. 8 1. The length of the trochlea, M4, is the length of the chord joining the crossing points of the anterior and the posterior ridges with the medial lengthwise curvature of the trochlea. 2. The width of the trochlea, M5, is the length of the chord P ipmr P iplr joining the medial and the lateral ridges. 3. The posterior width of the trochlea, M5(1), is the length of the chord P ppmr P pplr joining the medial and the lateral ridges. 4. The anterior width of the trochlea, M5(2), is the length of the chord P apmr P aplr joining the medial and the lateral ridges. 5. The medial chord of the trochlea, Cm, is the chord length between P ppmr P apmr . 6. The lateral chord of the trochlea, Cl, is the chord length between P pplr P aplr . 7. The medial arc of the trochlea, Am, is the length of the curve joining P ppmr and P apmr along the medial ridge. 8. The lateral arc of the trochlea, Al, is the length of the curve joining P pplr and P aplr along the lateral ridge.

These chord and arc measurements provide the support for the morphometric analysis performed at the reconstruction and the prototyping stage of talus repository reproduction.

On the right side in the figures, the standard deviation distribution for the sampled values is shown. It could be seen that range intervals of the standard deviations decrease for measurements on the 3D prints. In fact, evaluate original specimen with partially damaged ridges is subjective depending on the human operator interpretation. While, after the reconstruction, the corresponding parts on the boundary are "repaired" and thus the positioning of the characterictic vertices is less uncertain. Moreover, range intervals of the chord and arc lengths on the tali tend to increase for their counterparts on the 3D prints in average with 1.5mm. On the contrary, length estimations on numerical models oscillate around manually measured values within 1mm absolute value interval. One can conclude that the dominant of the error in the pipeline occurs during the reconstruction and consequently further improvement should be supplied.

Conclusion

The present work relates our experiences in teaching geometric modeling fundamentals through the laser scanner acquisition based pipeline and with a particular application domain of cultural heritage preservation. We choose to cover the pipeline as a whole starting with the original specimen, scanning, registrating and reconstructing the object boundary, and finally, printing it. The major advantage of this approach is that in addition to the classical visual feedback in the visual systems, a material feedback is also produced and in this way the correctness of treatments all over the pipeline is evaluated. The notion of error is illustrated in a very intuitive way through length measurements: manual ones with calipers and millimetre ruler band on the tali and their 3D prints, and numerical ones on the numerical reconstructions. All over the covered subjects, students are involved with challenging scientific problems looking for efficient computation solutions. Our believe is that this approach gives a deeper understanding of both theoretical and application issues in geometric modeling.

Fig. 1 .

 1 Fig. 1. Polygon Editing Tool: (left) Two scans (center) Scan registration (right) Scan integration

Fig. 2 .

 2 Fig. 2. Polygon Editing Tool: erroneous acquisition

Fig. 3 .

 3 Fig. 3. Left talus: (a) Superior view (b) Anatomical position (c) 3D print

Fig. 4 .

 4 Fig. 4. Anatomical landmark points on the trochlear surface.

Fig. 5 .

 5 Fig. 5. M. Fichoux and D. Sabary's programming project : geodesic path calculation.

Fig. 6 .

 6 Fig. 6. Ch. Mary and Ch. Delmotte's programming project : offsetting for 3D printing.

Fig. 7 .

 7 Fig. 7. J. Laviole and R. Leguay's programming project : discrete surface curvature evaluation.

  (a) for the original individual, in Fig. 8(b) for the reconstruction, and in Fig. 8(c) for the 3D print. By analogy, in Fig. 9(a), Fig. 9(b) and Fig. 9(c) the measurements on BCT92-S380-G are detailed. In each figure, on the left side, computed values are provided. The notations in use are as follows:

Fig. 8 .

 8 Fig. 8. Chord and arc measurements on BCT92-S319-G: (a) Original specimen (b) Numerical reconstruction (c) 3D print

Fig. 9 .

 9 Fig. 9. Chord and arc measurements on BCT92-S380-G: (a) Original specimen (b) Numerical reconstruction (c) 3D print