Mickaël Péchaud Di

Gabriel Peyré

Renaud Keriven

Extraction of Vessels Networks over an Orientation Domain

This paper presents a new method to extract a network of vessels centerlines from a medical image. The network is composed of local geodesics over a four-dimensional space that includes local orientation and scale. These shortest paths follow closely the center of tubular structures and can deal robustly with crossings over the image plane. The vessel network is grown by an iterative algorithm that distributes seed points according to a geodesic saliency field. Numerical experiments on a database of synthetic and medical images show the superiority of our approach with respect to several methods based on shortest paths extractions. With a minimum of user interaction, it allows to compute a complex network of vessels over noisy medical.

Introduction

Detection of vessel networks in images is of primary interest to help medical diagnostic. The extraction of an accurate network allows one to compute meaningful information such as the local width of the vessels and the connectivity of the networks from a single planar observation.

Crossings and occlusions in 2D medical images are challenging problems that cannot be overcome by a local search procedure. This paper recasts the 2D extraction problem as a geodesic computation over a 4-dimensional space. An additional scale dimension gives access to the local width of the vessels, and allows the direct extraction of the centerline of the vessel. A rotational dimension reduces erroneous detection when two vessels are overlapping. Numerical experiments show the superiority of our approach with respect to state-of-the-art geodesic based methods. A network extraction method relying on this framework is then proposed, and experimented on several real and synthetic medical images.

Previous Work

The problem of vessels extraction has received considerable attention in the computer vision and medical imaging communities -see for instance the recent review by Kirbas and Quek [24], or [16] for further references. Several class of methods have been proposed to segment vessels in medical images. They generally rely on the use of a local detector, postprocessed by a method that links locally detected vessels.

Local detectors allow to detect points belonging to vessels or portions of vessels depending on the modality of the image. Local detectors include various methods: thresholding of images intensities, ridge or crest detection [1,36,5], wavelets [21,46,10], differential operators [33,27], gabor filters [38] or matching filters [15, 7, 20, 3, 6, ?, ?].

Many methods allow to link or post-process the locally detected points. Among classical methods (inspired by ideas which early arose in computer vision community for edge detection [4]), thresholding [50], fusion processes [23,17,31,53], region growing algorithms [13,42,55,19], front propagation [30,44], or pixel classification [10,39,45] techniques were proposed.

Geodesic based methods is another class of methods allowing the linking of local featuresusually pixels intensities: the notion of shortest path or geodesic proved to be efficient for the extraction of salient curves in 2D or 3D images, see for instance [9,34]. Geodesic curves can also be used to extract tubular structure centerlines in 3D, as proposed by Deschamps and Cohen [12] and by Santamaria et al [40]. Li and Yezzi [28] proposed to extend the shortest path computation to a higher dimensional domain. They include local radius of the vessels as an additional scale dimension in order to stabilize the computations and to select the centerline without any post-processing.

Another way of linking local features is the class of tracking methods which start from a point belonging to a vessel (either user-defined, or detected using a ad-hoc method with respect to the modality), and iteratively track the vessel by analyzing the neighborhood of the current point in the direction of the vessel (look-ahead) [37,26,41,14,29,56,51,11,49]. While some of these methods can handle junctions, they usually fail to deal with crossings -which are an extra difficulty of 2D images.

Our method operates a mix of these categories. It lifts the 2D image in a 4D orientation space using local detectors of vessels at different orientations and scales -the 4D orientation space disambiguates crossing configurations [22]. The method is independent from the local detector used, which can be tuned precisely to the targeted application. Our method then uses a geodesic based formalism to compute optimal paths in this 4D space, leading to a robust global segmentation of vessels. Unlike methods which rely on a post-processing skeletization to compute the centerlines of the vessels [52,47,18], our method directly and naturally computes both centerlines and radii of vessels. Finally, we propose to use this segmentation to design an algorithm for network extraction. Based on a tracking on extended neighborhoods, our algorithm handles difficult crossing configurations.

[2] recently proposed a similar growing-of-minimal-paths framework, but it is specialized to the segmention of closed curves in 2D and meshing of surfaces in 3D.

Overview

Our vessel network is computed over a 4D domain that includes both local orientation and scale information. The metric over this lifted space is described in Section 2.

Vessels are computed as shortest paths over the 4D domain. The computation of these geodesics is exposed in Section 3. The performance of this method is evaluated in Section 4 on both synthetic and medical images.

We then define a network as a union of vessels. It is computed by an iterative procedure described in Section 5. Numerical examples on synthetic and medical images (both cortical optical imaging [25] and retinal imaging [32]) show the ability of our method to extract network.

Local Scale and Orientation Detection

Local Vessel Model

A medical image is treated as a 2D function I : [0, 1] 2 → |. The local geometry of a vessel is captured with a vessel model M (x) ∈ R for x = (x 1 , x 2 ) ∈ Λ = [-Λ 1 , Λ 1 ] × [-Λ 2 , Λ 2 ].
This model is a 2D pattern that incorporates our prior knowledge about both the cross section of the vessels and the regularity of vessel.

The prior on the cross section of the vessel is included by considering models M (x 1 , x 2 ) = m(x 2 ) that only depends on a 1D profile m. The prior on the regularity of the vessels corresponds to the ratio Λ 1 /Λ 2 of the horizontal and vertical dimensions of the model.

Model cross-section for vessels. A 1D profile adapted to both cortical optical imaging and retinal imaging is defined as

m(x 2 ) = 1 for |x 2 | > Λ 2 /2, exp(-α (1/2) 2 -(x 2 /Λ 2 ) 2 ) otherwise.
(1) The value α ≈ 0.05 was evaluated from a set of typical cortical images. Section 4.1 shows that our vessel extraction method is robust with respect to approximate choices of this absorption parameter.

Regularity selection. The ratio Λ 1 /Λ 2 of the model dimensions acts as a prior on the regularity of typical vessels. The more typical vessels are curved, the smaller Λ 1 /Λ 2 should be. Also, robustness to noisy images forces to use a model with a large enough area Λ 1 × Λ 2 . The value of (Λ 1 , Λ 2 ) = (1, 2) is used in our numerical experiments. It was selected to give the best vessel detection results on a set of exemplar images.

To overcome the inherent difficulties of the 2D detection problem, additional scale and orientation dimensions are introduced to increase the detectability of vessels.

Rotated and Scaled Models

The normalized pattern M (x) is rotated and scaled to match the varying orientation and width of vessels. Beside the choice of the pattern cross section m and the dimension Λ 1 × Λ 2 of the model M (x), the scaling of Λ(r) with r is another avenue to introduce some prior about vessels in the image. Small scales cortical and retinal vessels are less regular than large scale vessels. We thus chose to scale the dimensions of the model Λ(r) = rΛ linearly with the radius r. This causes thin tubular structures to be detected using a finer correlation analysis.

The warped model M r,θ (x) for x ∈ Λ(r, θ) = R θ Λ(r) is defined as

∀ x ∈ Λ(r, θ), M r,θ (x) def. = M (R -θ (x/r)) (2)
where R θ is the planar rotation of angle θ.

Figure 2 shows examples of models defined with (1) that are rotated and scaled according to (2).

Scale/Orientation Lifting

The image I is lifted in a 4D space by adding a scale and an orientation dimension. Let Ω be defined by

Ω def. = [0, 1] 2 × [R min , R max ] × [0, π) (3)
We call lifting the function F computed as the normalized cross-correlation between the image and the local model ( 2)

∀ ω = (x, r, θ) ∈ Ω, F (ω) def. = NCC Λ(r,θ) (M r,θ (•), I(x + •)) (4)
where I(x+•) is the image translated by x, NCC A (f, g) is the normalized cross-correlation between f and g over the domain A, defined by:

NCC A (f, g) def. = A (f -f )(g -ḡ) A (f -f ) 2 A (g -ḡ) 2
(5)

where h = ( A h)/|A|, |A| being the area of A. This lifting separates real 3D vessels that overlaps when projected at the same location by the imaging system but have different orientations.

The value F (x, r, θ) ranges from -1 to 1 and measures the likelihood of observing a vessel at a given location x with a width r and an orientation θ. The normalization of the detector makes it invariant under to intensity variations that occurs in medical images due to the elevation variation of the vessels and the imperfection of the imaging system. Adding a scale dimension yiels a robust and regularized estimation if the radius and the center of vessels.

Figure 3 shows an example of a cortical image where orientation lifting is crucial to distinguish locally between orientations. Numerical computations. A medical image is acquired on a on discrete grid of n × n pixels. The 4D lifting is computed for n r radii evenly spaced in [r min , r max ] and n θ orientations evenly spaced in [0, π), with n r = 12 and n θ = 12 in the experiments. This requires O((r max n) 2 n 2 n r n θ ) operations with r max ≪ 1 and n r , n θ ≪ n.

Geodesic 4D curves

Lifted Metric

The 4D lifting (4) defines an isotropic metric ρ over the 4D domain Ω

∀ ω ∈ Ω, ρ(ω) = max(1 -F (ω), ε). (6)
The parameter ε prevents the metric to vanish and is set to ε = 10 -3 in the numerical tests. This metric encodes local information about the presence of a vessel at a given position, scale and orientation. The extraction of minimal paths over this metric is a global optimization that links this local information and robustifies vessels extraction.

Distance Map and Geodesic Computation

The length of a lifted curve γ : [0, 1] → Ω over the lifted domain is defined as

L F (γ) def. = 1 0 ρ(γ(t))||γ ′ (t)||dt. ( 7
)
where the length of the speed vector

v = γ ′ (t) = (v x , v r , v θ ) is ||v|| 2 def. = v 2 x + λv 2 r + µv 2 θ , (8)
(λ, µ) being normalizing constants that controls the penalty on scale and orientation variations along the vessels in the images. They are set to λ = 0.5 and µ = 0.1 in numerical experiments. In practice we observed strong robustness with respect to the choice of (λ, µ). Given a set A ⊂ Ω of seeds points and a set B ⊂ Ω of ending points, a shortest lifted curve γ * (t) ⊂ Ω joining A to B is defined as a shortest path for the metric

L F γ * (A, B) = argmin γ∈π(A,B) L F (γ), (9)
where π(A, B) is the set of curves γ such that γ(0) ∈ B and γ(1) ∈ A. The corresponding geodesic distance is d F (A, B) = L F (γ * ). This definition can be specialized to a single starting point A = {ω 0 } and/or to single ending point B = {ω 1 } to define the geodesic distance between points and/or sets, e.g. d F (ω 0 , ω 1 )

def.

= d F ({ω 0 }, {ω 1 }). In practice, γ * is estimated as follows [44]: the distance to the seeds A is the geodesic action map U A (ω) = d F (A, ω) is the unique viscosity solution of the Eikonal equation

∀ ω ∈ Ω, ||∇U A (ω)|| = ρ(ω), with ∀ ω ∈ A, U A (ω) = 0 (10)
where

∇U A = ∂U A ∂x , λ ∂U A ∂θ , µ ∂U A ∂r T
, and where one has to be careful about computing the derivatives with respect to θ modulo π.

Calling ω 1 ∈ B the point in B with lowest distance to A, the geodesic curve γ * between A and B is then obtained by a gradient descent of

U A dγ * dt (t) = -∇U A (γ * (t)) with γ(0) = ω 1 (11)
Numerical computations. The numerical computation of U A on a discrete grid of N def.

= n 2 n r n θ 4D-points is done in O(N log(N )) operations with the Fast Marching algorithm of Sethian [43]. This algorithm solves an upwind finite difference discretization of (10) with an optimal ordering of the grid points. This ordering corresponds to a front propagation that starts from A and progressively propagates the values of U A in increasing order. The orientation dimension θ ∈ [0, π) is handled with periodic boundary conditions 0 ≃ π. A first order Euler-scheme was used to perform the gradient descent. Note that γ * is computed with sub-pixel accuracy.

Shortest Paths and 4D curves

A 4D curve c(x, x ′ ) between two points x, x ′ ∈ [0, 1] 2 is computed as a 4D geodesic in Ω between the 4D lifted sets A(x) and A(x ′ ) defined as

A(x ( ′ ) ) def. = (x ( ′ ) , r, θ) \ r ∈ [r min , r max ], θ ∈ [0; π) . (12)
The 4D curve is then defined

c x,x ′ def. = argmin γ∈π(ω,ω ′ ),ω ( ′ ) ∈A(x ( ′ ) ) L(γ). (13)
This 4D curve contains three components c x,x ′ (t) = (x(t), r(t), θ(t)). The path x(t) ⊂ [0, 1] 2 is the actual centerline over the image plane, whereas r(t) and θ(t) give the local width and orientation of the vessel, see Figure 4. These orientation and scale are important in medical applications, for instance to compute the local blood flow or to detect tissue diseases.

Evaluation of the Geodesic Centerlines

The accuracy of the centerline extraction is compared to two other methods on both synthetic and real data.

• In [8], a 2D metric is computed from an image intensity blurred with a Gaussian filter. The filtering helps to re-center the geodesic since the smoothed image exhibits a local maxima around the center of the vessels, at the cost of a loss of spatial resolution. The smoothing intensity is optimized to give the best results in each test.

• In [28], a 3D (space+scale) metric is computed. All the parameter of this 3D model are optimized to give the best results.

In the experiments, our method was used with an absorbtion parameter α = 0.1, which is not optimezed to fit the α of real or benchmark images.

Evaluation on Synthetic Crossings

The 4D lifting (4) is challenged by testing the extraction of a curved vessel with a self crossing. Figure 5 shows that the vessel curve is not correctly extracted with a metric that does not take into account the local orientation. A 2D purely spatial metric extracts a curve that is not well centered and does not capture the correct topology of the vessel. A 3D space+scale metric correctly centers the geodesic in the vessel but fails to detect the crossing singularity. Our 4D centerline position+scale+orientation favors the extraction of a longer curve that is both well centered and geometrically faithful to the true vessel.

Accuracy and robustness to noise

The precisions of the three algorithms are tested on several phantoms images, for which the centerlines positions and radii have known analytical forms with sub-pixelic accuracy. The cross section of these phantoms corresponds to the model (1) with parameter α = 0.01. An additive Gaussian white noise with various amplitudes are added to the phantoms. Ten phantoms are generated for each condition, and each noise level, see Figure 6. In each case, the true starting and ending points of each phantom are used, as well as the true starting and ending radii for the [28] method.

2D metric [8]

3D metric [28] our method: 4D metric (position) (position+scale) (position+scale+orientation) where t * is such that x * (t) is the ground truth centerline point closest to x(t).

Figure 7 shows Error C (c) and Error R (c) curves for several synthetic images as a function of the noise level. Using the 3D space+scale lifting [28] produces results of varying quality, and requires a careful tuning for each noise level of the parameter to achieve the optimal error rate. [8]with an optimal smoothing provides a precise evaluation of the centerline locations, but without any evaluation of the local radius. Our method provides both positions and radii with more robustness and accuracy.

Robustness to vessel model The dependency of our 4D lifting method with respect to the precise definition of the model M (x) is challenged by testing several values of the α parameter (1).

Table 1 reports the estimation errors Error C and Error R for phantoms f generated with α = 0.01 and a noise variance σ = 0.2. The centerlines extracted with models M (x) are computed for a parameter α set to three different values. This shows that the proposed method is able to recover with good accuracy the centerlines positions and radii, even if the vessel model is unknown.

Evaluation on Medical Images

Figure 2 :

 2 Figure 2: Vessel models for different orientations and scales. Here, Λ 1 /Λ 2 = 1/2 and m(.) is given by (1).This model encompasses medical knowledge about the light reflexion around blood vessels in cortical imaging. The image intensity inside a vessel is assumed to result from a light absorption (with coefficient α) proportional to the vessel width at this point. It is also widely used in the retina image community[7].The value α ≈ 0.05 was evaluated from a set of typical cortical images. Section 4.1 shows that our vessel extraction method is robust with respect to approximate choices of this absorption parameter.

Figure 3 :

 3 Figure 3: left: An original 2D image. right: Its 4D lifting (fixed radius), ranging from -1 (black) to 1 (white). White values indicate likely positions and orientations of vessels.

Figure 4 :

 4 Figure 4: Left: centerline extraction of a vessel in a cortical image. Starting point: white square. Ending point: black square. Right: corresponding orientation θ(t) and radius r(t).

Figure 5 :

 5 Figure 5: Comparison of the 2D [8], 3D [28] and 4D lifting (our method) when encountering a self-crossing.

Figure 6 :

 6 Figure 6: Some of the sample phantoms used in our benchmark (basic intensities range from 0 to 1), shown here with a spatially independent Gaussian noise of variance 0.15.The extracted 4D curve c(t) = (x(t), r(t), θ(t)) is compared to the ground trust c * using the following errors:Error C (c) 2 = 1 0 ||x(t) -x * (t * )|| 2 dt Error R (c) 2 = 1 0 |r(t) -r * (t * )| 2 dt (14)

Figure 8 ,

 8 Figure8, left, shows vessels extraction for a complex optical imaging of the cortex with several branches and intersections. As expected, the centerlines computed from different ending points are overlapping.


The granular Brownian ratchet: carrying out a historical thought experiment. We present an experimental realisation at macroscopic scale of the storied Brownian ratchet, which is an illustration of the Maxwell's demon. In our mechanism, the rotation of a centimetrescale 1D Brownian object in a granular gas is detected by an electromechanical converter (dynamo), generating a voltage proportional to its angular velocity. The current generated by this random rotation is rectified by an electronic device (demon), such that only positive current passes. Eventually, work is produced. The advantage of such a macroscopic setup is to allow measurement of all the observables per time unit: useful power (work), heat taken from the bath, and finally the efficiency of the equivalent heat engine. The feedback allowing the conversion from heat into work, expresses as a bias on the Brownian motion.

I. INTRODUCTION

The phenomenological theory of thermodynamics was founded in the XIX th century on primal principles, inspired by experience. Beyond the plainness of their statements, they reveal a profound significance [1]. The 2 nd law, among many other consequences, implies that heat spontaneously flows from high to low temperature bodies. James K. Maxwell was one of the first to sense the statistical nature of these axioms [2]. More precisely, he had the intuition that they should be taken into account only for the averaged observables. We now know that state variables are fluctuating quantities, but it was not the case then. Indeed, the conception of matter as composed of atoms and molecules in random thermal motion was not yet accepted in Maxwell's time. Such considerations were revolutionary and among the very first openings to statistical physics! To illustrate his point, Maxwell proposed an experiment that would, if performed, contradict the 2 nd principle, i.e., a paradox. The principle is as follows. If an imaginary being is able to take information on the state of a system (measurement), it can act on this system. For instance, it can sort particles or extract work from heat, at no cost of energy. This being has been later called Maxwell's demon. Different versions of this simplistic thought experiment fueled rich and profound discussions for over 150 years. The most famous are the Szilard engine [3], the Brillouin's rectifier [4], and the ratchet and pawl [START_REF] Feynman | The Feynman Lectures on Physics[END_REF]. The detailed operation of these devices cannot be discussed here, we refer to the detailed review from H. S. Leff and A. F. Rex [START_REF] Leff | [END_REF]. The paradox of Maxwell's demon has been first resolved by Smoluchowski, arguing that no sorting is possible if the demon itself is not at a temperature lower than that of the gas [START_REF] Smoluchowski | [END_REF]. In other words, the Maxwell's demon is actually a heat engine, the cold bath being not apparent at first sight (detailed in Feynman [START_REF] Feynman | The Feynman Lectures on Physics[END_REF]). Looking closer, the Maxwell's demon is operating a conversion from in- * antoine.naert@ens-lyon.fr formation into heat (entropy transferred to a cold bath). The detailed processed was first clarified by Brillouin [4], and Landauer [8], paving the way to the very active field of information thermodynamics [9][10][11]. Experiments were technically unrealisable in Maxwell's time. They could only be performed in the last decade, thanks to a colloidal Brownian rotor under feedback controlled rotating electric field [12], or a colloid bead oscillating randomly in the simple or double-well potential of optical tweezers [13]. These pioneering experiments made it possible to investigate quantitatively the functioning of Maxwell's demons. We present a novel experimental construction of the celebrated 'ratchet and pawl' [START_REF] Feynman | The Feynman Lectures on Physics[END_REF]. This is a mechanical realization of a Maxwell's demon, illustrating the arguments of Smoluchowski. In [START_REF] Feynman | The Feynman Lectures on Physics[END_REF], the random torque caused by the collisions with a molecular gas sets a rotor into a 1D Brownian motion. Consider a rectifying mechanical device, like a ratchet and pawl, fixed on its shaft, allowing the rotation in one direction only. In this case, work can be extracted, seemingly from a single heat reservoir, which constitutes a perpetual motion of the 2 nd kind. This is the Maxwell's demon paradox. We present here our first measurements in such system, of time-dependent observables such as heat and work per time unit. The originality of our setup is twofold. First, it is macroscopic (centimetre scale). The role of the molecular heat bath is played by a granular gas. Being dissipative, external power has to be supplied to sustain a Non-Equilibrium Steady State. It is a thermostat, in the sense that it sets an effective temperature [START_REF]Throughout this work, the temperature, kT , is given in energy units[END_REF], kT 1.36 10 -7 J [16]. Second, it makes use of an electromechanical conversion: the movement of the rotor is converted into an electrical voltage. Using an electronic rectifying device as the demon avoids at once all difficulties linked to mechanical plays and friction.

II. EXPERIMENT

A centimetre-size rotor is embedded in a granular gas playing the role of heat bath (about 340 stainless steel beads of 3 mm diameter and 0.1 mg mass). Such a macroscopic gas is dissipative, ie., some power is supplied to balance dissipation on the average, and achieve a steady-state. The beads are contained in a 5 cm cylindrical vessel, accelerated vertically at a few g (gravitational acceleration) by an electromagnetic shaker (B&K 4809), at a frequency of 40 Hz (see Fig. A rotating 20 × 20 × 0.4 mm stainless steel blade, immersed in the granular gas, is secured on the shaft of a DC micro-motor (Maxon RE10 118391), set vertically on the axis of symmetry of the 5 cm diameter cell. The micromotor is simply used as an electro-mechanical converter (dynamo). The voltage at its terminals is the image of the angular velocity: e(t) = α θ(t), thanks to Faraday induction. The prefactor α 0.0095 V/(rad/s) is characteristic of the device. The granular gas colliding the blade causes its random rotation, realizing a macroscopic 1D Brownian motion. If the electrical circuit is closed on a load, the fluctuating voltage at the terminals of the motor gives rise to a current i(t). This current is modulated by an electrical Maxwell's demon, composed of a voltage comparator and a switch. Whatever the voltage fluctuations, only one direction of current is allowed: the demon is an electrical rectifier, analog to the usual ratchet and pawl. To avoid spurious transients, we chose a high-speed voltage comparator, with compensated thermal drift (ref. Rt , where the total resistance in the loop is R t = R + r + r . The internal resistance of the motor is r 120 Ω. It is measured, like R, thanks to an ohmmeter. The current through the motor gives rise to a torque that opposes the rotation (Lenz's law). In other words, when the switch is closed, an Ohmic damping adds to the viscous damping of the gas. The ohmic load represents the weight lifted by the winch in the mechanical ratchet described in Feynman's lectures. The voltage noise e(t), represents the velocity of a rotor in the granular gas (hot bath). To detect properly the sign of the voltage, the thermal (Nyquist) noise of comparator itself must be much lower than e(t). This is obviously the case at room temperature (k B T 10 -21 J). This is, in electrical words, the expression of the Feynman-Smoluchowski's explanation of Maxwell's paradox: work can be produced because the heat engine is coupled to, both a hot bath and a cold bath. The demon stays in the low temperature heat bath.

III. THE DEMON'S CHOICE

Let us now compare the motion of the rotor with and without the demon. In the latter case, the rotor is subject to the random forcing of the gas, whereas the averaged collisions cause damping (gas viscosity). This equilibrium on average, is the essence of the fluctuationdissipation theorem, taking into account the temperature of the gas [16]. On average, the angle of the rotor remains constant (blue trajectory in Fig. 3). When the demon is at work, the symmetry between positive and negative rotation breaks, because ohmic dissipation in R occurs in one direction only. It results in a drift θ (temporal average). One example is presented in Fig. 3 Let us now consider the variance δ τ θ 2 -δ τ θ 2 of the angular time increments δ τ θ(t) = θ(t + τ ) -θ(t), over a time τ . It is the time coarse-grained mean square rotation. One clearly identifies In Fig. 4 a ballistic and a diffusive regime, where the variance is proportional to τ 2 or τ , for time-lags respectively smaller or larger than a certain crossover time. We note that this time is approximately the correlation time of the angular velocity: τ * 18 ms. Remarkably, δ τ θ 2 -δ τ θ 2 has exactly the same time dependance whether the demon is active or not: red and blue lines are exactly superimposed in Fig. 4. It means that the rotor undergoes a similar diffusion process, if considered in a frame of reference rotating at speed θ . Note that the average rotation with the demon is quite small, typically 5 rotation per hour, it is indeed a tiny perturbation onto the gas. Consistently, the rotor speed histograms are almost the same with or without demon (see Fig. 5). On the histograms, the negative fluctuations, when switch is open, are of course unaffected ( θ < 0). It could be different on the θ > 0 side, where the demon causes a very slight depletion, due to the damping of the blade motion by the load. The effect on the histogram of velocity is so small it can barely be seen. The overall effect is a slight angular drift under the action of the demon, already mentioned with Fig. 3.

The smallness of the perturbation, θ 0.016 rps, can be understood in different ways. It confirms that the fluctuations of the rotor are strongly coupled with the dynamics of the beads, and the energy taken from the bath is small. The bath is a good thermostat, for a demon that makes a weak perturbation. The measurements of hear and work rates are discussed in the next section.

IV. HEAT, WORK, AND EFFICIENCY

The paradox of the Maxwell's demon, as explained in The Feynman Lectures, is to consider the ratchet and pawl system as a heat engine, working between a hot bath on the one hand (here the granular gas at temperature kT ), and a cold bath on the other hand (the surroundings). The electronic demon system stands at room temperature (comparator and switch), as well as the load R and the motor. This model thermal machine extracts heat from the hot bath and releases part of it into the cold one, producing some work in between. Now, let us see the energy balance of this steady state heat engine. Per time unit, ẇ = Ri 2 represents the work produced thanks to the ratchet heat engine. On the other hand, the heat per time unit extracted from the hot bath is q = ei. Both histograms show fat tails for positive power fluctuations (see Fig. 6). FIG. 6. Normalized histograms of the heat and work per time unit, q (orange doted line) and ẇ (green line). Averaged heat and work per time unit are q = ei 6.6 µW and ẇ = R i 2 4.2 µW .

Two comments can be made: q and ẇ are positive, and exhibit a zero most probable value, corresponding to the switch opened. These statistics deserve specific investigation, but that is beyond the scope of this article. Let us consider the averaged quantities. Per time unit, the heat taken from the granular gas is q 6.6 µW, and the work produced is ẇ 4.2 µW. We note again that these values are very small if compared to the power used to keep the granular gas in a steady state, ∼ 10 W, called house keeping heat [17]. It is a reliable thermostat. The efficiency of the heat engine is η = ẇ q 64%. An upper bound for the efficiency corresponds theoretically to a reversible exchange. It is the Carnot efficiency: η C = 1 -kTc kT h , where kT c and kT h are the temperatures of the cold and hot baths, respectively. The temperature of the gas is kT h = 1.36±0.15 10 -7 J, and the room temperature is kT c = k B T room 10 -21 J. This leads to η C 1 -10 -14

1. As expected, η < η C . The difference is mainly due to the losses. To be specific, it is mostly Joule dissipation: (r + r ) i 2 2.3 10 -6 W. One must mention another contribution, the heat released to the cold bath because of the entropy production in the erasure of the demon's information, as proposed by Brillouin [4]. It can be obtained from the digital switching signal of the demon, k B T log 2 per bit of information erased. As it is, at room temperature, this contribution is extremely small, barely measurable (∼ 3 10 -18 W on the average). The investigation of this digital signal is in progress.

V. DISCUSSION

We designed and produced a macroscopic version of the sometimes-called 'Feynman-Smoluchowski ratchet', a thought experiment that allows to convert heat into work through the action of an information device, a Maxwell's demon. Our setup is operating in contact with a granular gas as heat bath. It makes use of, first, electromechanical conversion, and, second, an electronic device as demon, to replace the traditionally represented mechanical ratchet and pawl system. Granular ratchets have been studied in the past, but, to our knowledge, no measurements were performed on their thermodynamic properties [18]. Our system is profitable in that no friction exists in the demon's part. The only relevant damping effects are the finite impedance of the motor, and in the open and closed states of the switch, easy to take into account. The dissipation in the gas itself is evidently irrelevant. Indeed, as it is merely a source of noise, the dynamics of the rotor is not affected by the dissipation in the beads' collisions. Only the effective temperature matters. We deliberately chose Smoluchowski and Feynman's take on Maxwell's demon paradox. In that regard, it is a heat engine running between a hot heat bath (granular gas) and a cold one (room). We measured time series of the heat extracted from the hot reservoir, and the work performed on a load (here represented by R). Thanks to this macroscopic approach we can measure precisely various fluctuating observables, such as heat and work per time unit: q(t) and ẇ(t). We can therefore access the stochastic efficiency, defined as the ratio of the instantaneous work produced over heat consumed at the same time: η(t) = ẇ(t) q(t) . Such quantity, fluctuating in small scale systems, has attracted increasing interest in recent years [19]. One must note an experimental contribution on colloidal cyclic heat engines [20]. More interesting, our Brownian ratchet gives a direct access to the measurement of the entropy of information. It is a convenient experimental test rig, allowing to measure any fluctuating observable in a reliable way. Besides, from a more practical viewpoint, macroscopic Brownian ratchets typify the very principle of energy harvesting, ie., extracting useful work from fluctuations [21].

M

  . Lagoin, C. Crauste-Thibierge, and A. Naert * ENSL, CNRS, Laboratoire de physique, F-69342 Lyon, France. (depuis le 23/3/22) (Dated: April 24, 2022)

  FIG. 1. Photograph of the experimental system. The beads are vibrated vertically in a vessel by an electromagnetic shaker. A blade is immersed in the grains, fastened on the shaft of a small DC-motor. (Exposure time: 300 µs.)

FIG. 2 .

 2 FIG.2. Electrical circuit associated with the probe -The motor M is represented in a dashed rectangle by an ideal voltage source e and its internal resistance r. The demon D is composed of a comparator, at the terminals of the motor, triggering a switch that allows the current to pass or not in the load resistance R, depending on the polarity of u1.

9 FIG. 4 .FIG. 5 .

 945 FIG.4. Variance of the angular increments δτ θ over the timelag τ with demon (red line with symbols) or without (blue line). As an eye guide, doted lines of slopes 1 and 2 are drawn. The insert shows 1 τ ( δτ θ 2 -δτ θ 2 ) against τ . Vertical dashed lines figure the velocity correlation time, τ * 18ms.

ACKNOWLEDGEMENT

We acknowledge F. Ropars for designing and realising the electronic demon, and also J. Marot. We are grateful to S. Aumaître, A. Archambault, A. Steinberger, for many fruitful discussions, and V. Moncorgé for help in taking the photograph (Fig. 1). We also thank M. Adda-Bedia for careful reading and suggestions on the manuscript.


Yann Savoye

Alexandre Meyer
email: alexandre.meyer@liris.cnrs.fr

Workshop in Virtual Reality Interactions and Physical Simulation VRIPHYS (2008) F. Faure, M. Teschner (Editors

Keywords: (according to ACM CCS): I, 3, 7 [Three-Dimensional Graphics and Realism]: Animation, Virtual reality; I, 3, 5 [Computer Graphics]: Computational Geometry and Object Modeling

Real-time animation of human-like characters has been an active research area in computer graphics. Nowadays, more and more applications need to render various realistic scenes with human motion in crowds for interactive virtual environments. Animation and level of detail are well explored fields but little has been done to generate level of detail automatically for dynamic articulated meshes. Our approach is based on the combination of three interesting layers for run-time level of detail in character crowd animation: the skeleton, the mesh and the motion. We build a Multiresolution Skeletal Graph to simplify the skeleton topology progressively. In contrast with previous works, we use a Dual-Graph Based Simplification for articulated meshes, where the triangle decimation is driven by triangle compactness, to build a dynamic, continuous, progressive and selective mesh level of detail. We also present Power Skinning to ensure the stability of Linear Smooth Skinning, during the simplification, with an efficient multi-weight update rule.

Introduction

Real-time animation of human-like characters has been an active area of research in computer graphics for several years. In addition to crowd behavior simulation, which is not the topic of this paper, the independent animation of each character inside of a large crowd in real-time poses a major problem in terms of performance and realism. Motion capture data sets and animation techniques are used widely in a variety of applications, such as medicine, sports, serious games, video games, and film production, where synthetic crowd generation and processing techniques are needed. In this paper, we present an overview of the related works in a number of different areas in order to identify the weaknesses of these methods and to present our contributions. We introduce a new three-layer level of detail hierarchy model specific to bone-controlled articulated meshes. A critical analysis of work performed and results achieved allow us to think about new perspectives for future work. This work has been realised in the context of the PlayAll project.

Computing the mesh deformation of every character is CPU time consuming. Moreover, avoiding repetitive animation is hard to do by hand for an artist. The principle objective of this paper is to propose an automatic method for generating relevant level for animated characters by adapting, at run-time, the complexity of the models without resorting to pre-computation. Indeed, it is not necessary to display or animate all details of the geometry of a character when it is far from the observer. Real-time display is possible if the character animations and models are set at different levels of detail: a precise level for characters close to the camera and more coarse levels for more distant characters. Incontestably, the generation of animation on a large scale needs to take level of detail into account inside the graphics pipeline.

Background and Related Works

The basic idea is to view character animation as an articulated mesh or skeleton-driven animation problem. We propose an overview of fundamental tools used in character animation and level of detail that are useful for crowd animation.

Skeletal Animation

The realistic rendering of human models uses a simple anatomical character layered model with four layers: skeleton layer, muscle layer, skin layer and cloth layer. In our work, we define the character model using a deformable, articulated layered model with two layers: skeleton and skin. Most real-time characters are modeled as an arbitrary 2manifold triangle mesh that provides volume to the skeleton and simulates the skin. Every mesh is defined by its geometry and its connectivity. A skeleton of animation S = (J , F) is made of a set J of several joints assembled in a hierarchy and a set F of frames called the motion data. A skeleton may be thought of as a collection of rotational joints. Motion data consists of a bundle of motion signals sampled in a sequence of frames.

The mesh M can be attached or bound to an underlying skeleton through a process called skin binding. Each vertex of a polygonal mesh is attached to one joint or more joints. A blend weight w j k ,vi can be defined to measure the influence of a joint j k on a vertex v i . We assume that w j k ,vi = 0 if the vertex v i is not influenced by the joint j k . We define the complete joint-vertex influence set by a list of tuples to represent the joint-vertex attachments or assignments:

A = j k ; v i ; w j k ,vi , ∀ j k ∈ J , ∀v i ∈ M
The full joint-vertex influence could be estimated and adjusted by computing normalized vertex weight values for each association between a vertex and a joint. The Point-toline Mapping, Containment Binding and Nearest Joint Assignment algorithms can be used to adjust weights when each vertex is attached to a joint. Linear Smooth Skinning (LSB) is a widely used and important traditional skinning algorithm known under many names, such as Smooth Skinning and Skeleton Subspace Deformation (SSD). This technique works by attaching a mesh vertex v i with a set of joints. The sum of all weights for a given vertex v i must be normalized.

n ∑ k=1 w j k ,vi = 1 ; 0 ≤ w j k ,vi ≤ 1
A blend vertex is obtained by determining the result transformation of the weight average of all attached joint transformations. We obtain the follow formula:

v i, f global = n ∑ k=1 w j k ,vi • (v i ) local j k , f •W j k , f
where:

• (v i ) local j k , f
is the default position of the vertex v i in the j k joint local coordinate system deformed by the frame f ; • W j k , f represents the world space matrix associated with the joint j k for the frame f ; • w j k ,vi represents the normalized weight for the joint j k associated with the vertex v i .

Level of Detail Overview

The computation and storage requirements for scenes such as those used in virtual reality applications by far exceeds the capacity of modern graphical hardware. a motion level of detail framework.

Weaknesses in current methods

Weaknesses in current methods often appear in video games that use a reduced collection of pre-computed static meshes for the same character at various levels of detail. LOD techniques based on mesh simplification like [START_REF] Hoppe | Progressive meshes[END_REF] can produce low detail approximations for a single pose but do not provide any guarantee of fidelity for other poses leading to animated characters that appear non-continuous and unrealistic. Techniques, like the Quadric Error Metric, that work well on a single character, are effectively impossible to apply in the context of a crowd scene where characters are designed to be unique. Ensuring the best quality of a simplification requires a huge computational cost for no significant visual improvement when a character is far away, masked partially by other characters or elements of the decor.

Contributions

We introduce a multi-layer LoD hierarchy model specific to articulated meshes by identifying three layers for level of detail in character animation. Our approach called Multi-layer Animation Level of Detail adaptively changes the accuracy level of each character in a crowd in real time, handling bottlenecks through several primary contributions including:

• Skeleton Labeling Algorithm: we suggest a new algorithm to label joints and to cluster them by their importance in the topology in order to be able to sort them in a decimation order according to the joint energy. • Multiresolution Skeletal Graph: we present a dynamic skeleton simplification algorithm, based on a statistical estimation of the energy, which respects the homotopy of the original skeletal structure. A hierarchical skeleton at different resolutions is generated for the multiresolution skeletal graph. The correctness of the result is guaranteed by graph theory. • Dual-Graph Based Simplification for Articulated-Meshes: in contrast with previous works, our dynamic skeletal-controlled simplification of the mesh uses a new structure called BindGraphMesh embedding the dual graph of the mesh to pilot its simplification and embedding the vertex-joint influence set to update the binding information at simplification run-time with a multiweight update rule. In addition, we implement a simple, selective, progressive, dynamic, view-dependant and pose-dependent mesh simplification method to enhance the performance. • Power Skinning Algorithm: we present a new advanced method called Power Skinning for skeleton-driven control of the level of detail based on a multi-weight update rule for Linear Smooth Skinning to ensure and maintain the skinning performance during the dynamic mesh simplification.

Multi-layer Animation Level of Detail

For large virtual environments, such as a virtual city populated by human-like characters, only a relatively small proportion of characters are relevant to a viewer's interactions at any given moment. The rendering bottleneck of processing a highly detail geometric model is well-known in the computer graphics field. We identify three important layers of level of detail in a skeletal-driven character with skinned mesh vertices, where we can perform an automatic, selective, progressive, dynamic, view-dependent and pose-dependent simplification: the skeleton, the mesh and the motion. The first layer simplifies the skeletal animation by removing joints based on a joint energy metric. The second layer simplifies the mesh geometry by collapsing triangles. The third layer merges key frames to simplify animations. Finally, a Power Skinning algorithm is employed to ensure skinning performance during animation is stable.

c The Eurographics Association 2008.

Dynamic Skeletal Level of Detail

The idea of simplifying the skeleton is the very foundation of the use of skeletal animation. It was therefore decided to combine some bones and articulation, to get a simplified version that keeps enough degrees of freedom to represent the main movements. Skeletal simplification plays an important role in crowd animation; by simplifying the skeleton, we can reduce the cost of transforming joints. Previous research on skeletal simplification has focused on the improvement of performance of physically based simulation. We propose to exploit joint energy to control the decimation of joints for the progressive simplification of the skeleton. The principle of our method is to isolate sub-kinematics chains (for instance legs) by identifying key joints in the skeleton topology to decimate progressively intermediate joints according to the joint energy. We propose an overview of our Dynamic Skeletal Simplification Algorithm as follows:

• Step 1: loading the animation skeleton and calculation of all relative and absolute joints matrices (Figure 4

Skeleton Labeling Algorithm

The importance of a joint in the topology can be evaluated by identifying its type. A joint cluster represents a specific kind of joint according to its impact and its role in the topology and the motion. The Skeleton Labeling Algorithm is based on joint statistical study of skeleton topology with a computation of average values. We can define many categories in order to represent groups of joints with the same topological properties. In particular:

• χ: joints belonging to the end-effector joints group,

• ζ: joints belonging to the root joints group,

• ψ: joints belonging to the major joints group,

• κ: joints belonging to the minor joints group (unique removable joint group), • µ: joints belonging to the extreme joints group.

The joint statistical measured values for all joints j k in a given skeleton topology J is based on the information and statistical criterion that we detail here:

• ρ ( j k ) is the geodesic distance between j k to the root joint.

It corresponds to the depth of the joint in the tree hierarchy of joints, • neighbor ( j k ) is the total number of joints directly linked with it, • E ( j k , m i ) is the joint energy according to a given motion that we will discuss in details in the next subsection.

We also define functions to compute the average for all statistical criteria values. For example we define a function averagengb (J ) which takes a skeleton topology and computes the average of joints in the direct neighborhood for all joints and a function averagerootdist (J ) which takes a skeleton topology and computes the average geodesic distance between joints and the root joint:

averagengb (J ) = ∑ ∀ j k ∈J neighbor ( j k ) |J | averagerootdist (J ) = ∑ ∀ j k ∈J ρ ( j k ) |J |
The Skeleton Labeling Algorithm is defined as a heuristic based-rule mechanism on joint properties. This algorithm allows us to identify kinematics sub-chains by finding their extremity represented by two non-minor joints. For all joints j k ∈ J , we compute a label L ( j k ) ∈ {κ, µ, ψ, χ, ζ} using the following identification sort-rules:

1. If (ρ ( j k ) = 0) ⇒ (L ( j k ) = ζ) 2. Else if (childrennum ( j k ) = 0) ⇒ (L ( j k ) = χ) 3. Else if (neighbor ( j k ) > 1.9 * averagengb(J )) and (ρ ( j k ) < averagerootdist (J )) ⇒ (L ( j k ) = ψ) 4. Else if (neighbor ( j k ) > 1.9 * averagengb(J )) and (ρ ( j k ) > averagerootdist (J )) ⇒ (L ( j k ) = µ) 5. Else L ( j k ) = κ c

Joint Energy

The importance of a joint in a motion is dependent on its level in the hierarchy and can be evaluated by computing a metric we call the joint energy. Moreover the rotation of a joint belonging to the minor group affects the motion very little. We must have a good grasp of the data structure that represents the structure of the human body in order to classify joints for a given motion. We use the following notation:

• h: the depth of the tree representation of the hierarchy of joints, • d ( j k ) the depth of the joint j k ∈ J the hierarchy of joints.

We define the weight factor of the joint j k ∈ J by:

ρ ( j k ) = h-d( j k ) h
A joint closer to the root will have a bigger impact on the rotation of its children than a joint closer to the leaves. We define the energy of a joint for a given motion m i with f frames. Each joint j k has three Euler angles r k (t) = x j k (t) , y j k (t) , z j k (t) which represent the relative orientations of the joint j k to its parent joint. We estimate the joint energy E ( j k ) of a joint j k for a given motion m i as follows:

E ( j k , m i ) = ρ ( j k ) • 3 ∑ a=1 f -1 ∑ b=0 R a j k (b + 1) -R a j k (b) Max R a j k
where R 1 j (t), R 2 j (t), R 3 j (t) are the three Euler angles of joint j k at the frame t. Max R a j k is the maximum Euler angle of the joint j k that can be rotated about the a th axis.

Multiresolution Skeletal Graph

We propose to build a Multiresolution Skeletal Graph to produce a level of detail articulation (Figure 4,d). We build a continuous skeletal level of detail method with heuristics and rules. We build a level of articulation by studying the joint energy that represents the importance of each joint, reflecting how much the animation would be damaged by removing it. We sort joints such that joints with higher importance come first. The Skeleton Labeling Algorithm is key to allowing us to identify kinematics sub-chains inside the skeletal graph. A kinematics sub-chain corresponds exactly to the shortest joint chain delimited on each extremity by a nonminor joint.

Multiresolution Skeletal Graph Construction: Skeleton topology is a graph. We convert the tree of bone to a multiresolution graph to produce the skeletal simplification with the conservation of skeleton homotopy. The simplification is based on the Skeleton Labeling Algorithm to preserve a homotopic simplification. The Multiresolution Skeletal Graph is constructed as a collection of multiresolution kinematics sub-chains. To ensure a topologically valid method, we build a transitive closure on each sub-graph of each kinematic subchain.

Joint Decimation: In the next step, joint decimation implies a skeleton simplification. We allow simplification only inside kinematic sub-chains. Each joint, considered as a graph node in the multiresolution structure, could be enabled or disabled, except for joints belonging to the end-effector group, the root group and the major joint group which are always enabled. In the multiresolution graph, joints that do not meet an energy threshold are disabled as well as edges having at least one disabled joint.

Skeletal Topology Simplification and Extraction:

We are able to extract a skeleton tree for a given LoD threshold by a graph pruning process based on the joint decimation and shortest edges in the graph. The extraction of a LoD skeleton at a given resolution is based on the extraction of a minimal spanning tree with the minimal edge cost in the graph. We notice that only minor joints could be decimated in the topology during the skeletal simplification to preserve important hub joint in the topology.

Skinned Mesh Level of Detail

A common polygonal mesh modeling is a shared list of vertices and a list of faces storing pointers for its vertices. This representation using an Indexed Data Structure is both convenient and efficient for many purposes, however in some domains it proves ineffective. A mesh simplification algorithm is presented to reduce the number of vertices in a dense mesh of triangles during the animation. Simplifying deforming meshes using a standard algorithm is error prone. We present a run-time mesh simplification algorithm to reduce the number of vertices in a dense mesh of triangles. We propose an overview of our Skinned Mesh Level of Detail as follows:

• Step 1: a highly detailed arbitrary triangulated mesh is loaded inside the BindGraphMesh structure, the embedded dual graph of the mesh is built and embedded in the BGM structure; the vertex-joint influences of each vertex are also encapsulated in nodes. The BGM structure elements (faces, nodes) are automatically enriched by relationship information (Figure 8,a). • Step 2: nodes which are near a joint are marked as ineligible when the corresponding influence is measured by a strong blend weight above a given threshold (Figure 8,b). • Step 3: at the same time, we compute the triangle compactness or fairness and build a sorted sequence of original faces, according to triangle compactness, which are eligible for the collapse (Figure 8,c). • Step 4: given a LoD threshold, we do a triangle-collapse on triangles that meet the threshold in the sequence of collapsible triangles, and we undo the trianglecollapse operation on others. The decimation process re-meshes and re-textures identified Mesh Patches. The sequence of collapsible triangles is enriched and updated with the appearance of new triangles created during the simplification iterations (Figure 8,d).

• Step 5: Power Skinning is applied during the re-meshing to ensure skinning stability by calculating blend weights with the update multi-weight rule for the new single vertex of integration according to the sets of vertex-joint influences of the decimated triangle.

Dual-Graph Based Simplification for Articulated-Meshes

We present a new structure called BindGraphMesh (BGM) corresponding to an augmented mesh structure with an embedded mesh dual graph and joint-vertex influences in the mesh nodes. The dual graph of a mesh is a key point for a connectivity-oriented mesh simplification. A Bind-GraphMesh consists of a dynamic multiresolution mesh built at run-time and made of an original base mesh at full resolution and a set of collapsible triangles at a given character configuration.

Skeletal-Controlled Mesh Simplification

We introduce a simplification method driven by the skeleton joints. We build an automatic method to find specific parts of the model. We identify ineligible nodes corresponding with areas away from the joints that have been drastically simplified. We mark as ineligible nodes according to a threshold based on the higher subset skin weight of each joint. Our technique allows a skeleton-driven selective LoD and a controllable decimation during the simplification where ineligible vertices are conserved in the mesh. Areas away from the joints will be dramatically simplified, nevertheless we preserve high-detail regions around joints which are significantly deformed during the simplification. For a given character C = (S, M, A), we denote by v i ∈ M a vertex of the articulated mesh of the character C. The blend weight average for a given articulation j k ∈ J of C is defined by:

η (α) = ∑ w jα ,v i |{wj α,vi }|
We build the set of eligible vertices for collapsing operations, denoted by Ṽ, as follows:

∀v i ∈ M, v i ∈ Ṽ ⇐⇒ w j k ,vi < (τ • η (k) , ∀ j k ∈ J )
where τ is a parameterizable coefficient. We use τ = 1 by default in our implementation. If v i / ∈ Ṽ, v i cannot disappear during the mesh simplification.

Triangle decimation

The triangle decimation process is based on the trianglecollapse algorithm. Removing a triangle implies removing all triangles sharing an edge with this triangle, so a total of four triangles are removed. To preserve mesh topology, it is better to decimate triangles with small angles first. The decimation technique iteratively removes faces from the mesh, re-triangulating the resulting hole after each step. During a triangle-collapse, all vertices of the triangle merge into a new vertex that is carefully positioned to ensure the simplification. To achieve a good result for low computation cost, we collapse triangles at their barycenter. Triangle compactness allows us to sort triangles for decimation and to ensure the correctness of the simplification for a low computational cost and good enough results compared to the Quadric Error Metric. The triangle compactness measure, Compactness f ace , is equivalent to the aspect ratio and is given by:

Compactness f ace = 4 √ 3a l 2 1 +l 2 2 +l 2 3
where a is the area of the triangle f ace and l i are the lengths of its edges. Taking into account the aspect ratio of triangles helps to select operations that improve mesh quality. A flat triangle has a compactness near zero and should be decimated in priority (Figure 9). To generalize the triangle decimation to n vertices, we define a set of eligible vertices to collapse by a tuple (v 0 , . . . , v i , . . . , vn) ⊂ Ṽ. We denote by v the new single vertex resulting of the aggregation of (v 0 , . . . , v i , . . . , vn), and we describe this polygon collapse operation by the contraction as follows:

P = v i , ∀i ∈ [0, n] , v i ∈ Ṽ → v
In the triangle decimation case, three vertices (va, v b , vc) are contracted in a new representative vertex v that must be chosen. This collapse is driven by contracting two edges in the dual graph of the mesh. To ensure maximum of quality for the mesh simplification, we choose the triangle barycenter to be the mesh collapse vertex v.

For each decimated triangle, we identify its associated Mesh Patch (Figure 12), or a sub-part of the mesh, composed of this triangle and all neighboring triangles sharing a node of an edge in the dual graph. A triangle decimation involves the deactivation of all triangles inside the associated patch and the re-meshing of the oldest neighbors (orange triangle) in the interior of the patch with new triangles. We associate its list of old triangles in the patch before the simplification and a list of new triangles in the patch after the simplification with a new contracted position.

Power Skinning Algorithm

One problem remains: the skin. The original skin refers to all the bones. As we do LoD, we would like to have a reduced skin as well at run-time with aggregation of blend weights. Articulated meshes use weighted vertices in order to smooth blending at the joint, so we have to ensure smooth blending during the re-meshing process. The Power Skinning Algorithm is a Multiple Weight Update Rule that ensures Linear Smooth Skinning stability during re-meshing. During the triangle collapsing and re-meshing phase, a node integration is done and it is necessary to update the weight and texture coordinate information at the same time.

D = n ∑ i=0 d i = n ∑ i=0 d euclidean (v; v i )
The impact of a vertex of the mesh belonging to the triangle to decimate is inversely proportional to the distance, so we compute an Inverse Distance for each vertex vα ∈ P, α ∈ [0, n]:

dinvα = n ∑ i=0 d euclidean (v; v i ) -dα
For the whole configuration, we compute a single Inverse Power Distance (Dinv) as the sum of inverse distance between v and vα ∈ P, α ∈ [0, n]:

Dinv = n ∑ i=0 (D -d i ) = (n -1) D
We introduce the concept of a Power Weight, denoted by W, computed for each vertex vα ∈ P, α ∈ [0, n] in order to control the importance of all vertex-joint influence tuples for a given vertex participating in the aggregation of weights:

Wα = dinvα Dinv = D-dα (n-1)D
We use an accumulator to compute the set of vertex-joint influences for the new vertex v. The concept is that when a triangle is collapsed in a new vertex, the weights of each of the vertices of the collapsed triangle are transferred to the new vertex in parts inversely proportional to their distance from it. The consistency of v's blend weights is kept using a function of the inverse distance weighting (Figure 13.d) according to the vertices participating to the aggregation.

v i → j k ; w j k ,vi , w j k ,vi ∈ [0; 1] v → v; j k ; n ∑ l=0 W l * w j k ,v l , ∀ j k ∈ J       w v, j1 • • • w v, j k • • • w v, jn       =              n ∑ l=0 W l * w j1,v l • • • n ∑ l=0 W l * w j k ,v l • • • n ∑ l=0 W l * w jn,v l             
The stability of multi-weights skinning is ensured because n ∑ k=0 w v, j k = 1. Now, the new vertex v can be deformed correctly through Linear Smooth Skinning. We demonstrate the stability of the skinning of a multi-weighted articulated figure. Whereas our purpose is real-time simplification, Power Skinning may also be used by artists to automatically generate LoD levels using our progressive mesh based geometry simplification support.

Motion Level of Detail

After skeletal and then mesh simplification, we propose the third layer of our complete LoD scheme: motion simplification. An animation is a sequence of key positions. The motion structure yields two relationships among sampled values, inter-frame and intra-frame relationships. If each frame is considered independently, then there could be an undesirable jerkiness between consecutive frames. Therefore, we take into account the inter-frame relationship as well. A character shown close to the viewer can be very detailed, but c The Eurographics Association 2008.

once it moves away from the camera, we cannot see the individual fingers of a character anymore. Our motion LOD system will stop animating the fingers automatically, but there can still be a lot of overhead costs from the finger nodes; they still have to be transformed. Also nodes that add small detailed meshes to the character still exist in the skeleton and waste processor power. We merge nearby postures to simplify the number of frames. We build a multiresolution motion data tree where two neighboring frames can be fusionned if the distance of similarity between them is below a given threshold. We achieve an extrapolation by removing a frame from among the two or by merging the two into a single frame, preserving the duration of frames on the result. The distance of similarity of two frames for a given skeleton can be measured by the energy needed to move all joints from one frame to the other. In Figure 14, line 0 is the initial animation with n key positions.

Controlling the Level of Detail

Our system proposes level of detail control in interactive animation systems to pilot the progressive simplification by eliminating unperceivable details of skeleton, mesh and motion data. The LOD manager picks the best LoD of skeleton, mesh and motion according to the visual importance. Polygons, motion frames and joints removed according to the distance between the character and the camera. If the character is currently at a distance d from the camera, we denote the current interpolation value as α(t); we define 3 interpolation values: one for the skeleton (t skeleton ), another one for the motion (t motion ) and another one for the mesh (t mesh ). We denote by ts and te the start and the end value of the interpolation respectively. Such an interpolation is easy to use, computationally cheap, guarantees C 0 continuity, and ensures that our system dynamically produces a continuous level of detail.

Results & Discussion

We implemented a prototype version of a real-time environment with our multi-layer level of detail. Tests of our implementation were performed and were obtained using an Intel Core Quad CPU workstation with a 2.0 Ghz processor, 2 Go of RAM and a Nvidia Geforce 8800GTS video card. The real time rendering is used to measure the performance of structure construction and run-time processes. The computational time needed to generate the Multiresolution Skeletal Graph for one character is less than a millisecond for all tests. Table 1 gives a performance summary of the overall run-time simulation. The computational time needed to create the BindGraphMesh structure (t initBGM ) and lower mesh level of detail through a series of 10 trianglecollapse operations (t 10 * decim ) is provided, according to the number of joints ( joints) and the number of triangles (triangles), for one character. Figure 15 shows the gain of using level of detail on a skeleton used in a crowd made of one to 250 skeletons structures, comparing the frame rate at full resolution and the frame rate obtained with a decimation of 39% of the total number of joints. Notice that for 50 skeletons, our system is two times faster than without LoD. As seen in Figure 16, by eliminating 64% of the triangles of the mesh of a character with 3928 triangles at full resolution, we can improve the frame rate by a factor of 1.4. Moreover, errors caused by the simplification of the mesh have no perceptual impact when the character is relatively distant from the camera. The multiresolution skeletal representation allows us to produce a good skeleton LoD and huge crowds of characters that are extremely detailed. By progressively simplifying the skeleton we can reduce the cost of transforming joints automatically. Our Dual-Graph Based Simplification for Articulated-Meshes, using an iterative triangle decimation allows us to compute simplified meshes for a cheaper cost. Desired real-time continuous level of detail can be dynamically extracted at run-time.

joints

Conclusion & Future work

Real-time animation of deformable objects is always a tradeoff between visual fidelity and computational complexity. A simplification method for crowd scenes is a more challenging issue. In this paper, we present various layers for the application of level of detail specifically designed to work with bone-controlled articulated meshes and their synergistic usage to enhance the overall performance as shown by the results. We describe how to exhibit sensitive key layers for the level of detail through useful associated tools for real-time interactive crowd environments.

Because of a significant amount of programming required, there are some limitations to our prototype system that could be addressed in future works. We select several obvious areas for future works. Swapping among skinning methods could be one other method for performing a skinning level of detail. Indeed, the number of bone influences per vertex may be a fourth layer to our system. Additional work needs to be done to perform a continuous Skeletal LoD and Mesh LoD interlaced. Moreover adding detailed motion clips, with hands and fingers, to an existing skeleton and mesh can increase the realism and should better show the capacity of our three layer LoD technique.

  weight

Figure 2 :

 2 Figure 2: Binding.

Figure 3 :

 3 Figure 3: Linear Smooth Skinning.

Figure 4 :

 4 Figure 4: Multi-resolution Skeleton simplification.

  , a). • Step 2: computing skeleton statistics including the joint energy (Figure 4, b). • Step 3: application of the Skeleton Labeling Algorithm to cluster joints into categories (Figure 4, b). • Step 4: building the Multiresolution Skeletal Graph. A transitive closure is computed on each kinematic subgraph of each kinematic sub-chain and all edges are weighted proportionally to their length (Figure 4, c). • Step 5: extracting a LoD skeleton from the Multiresolution Skeletal Graph according to a certain energy threshold (Figure 4, d).

  The Eurographics Association 2008.

Figure 5 :

 5 Figure 5: Joint Labeling.

  energy

Figure 6 :

 6 Figure 6: Joint Energy.

cFigure 7 :

 7 Figure 7: LoD Skeleton extracted.

Figure 8 :

 8 Figure 8: Multi-resolution mesh.

Figure 9 :

 9 Figure 9: Dual-Graph Simplification (left) and Triangle Compactness (right).

Figure 10 :

 10 Figure 10: Ineligible Vertices (in red).

Figure 11 :

 11 Figure 11: Decimation Sensitive Deformation.

Figure 12 :

 12 Figure 12: Dual-Graph Based Triangle Decimation inside of a re-meshing patch.

Figure 13 :

 13 Figure 13: Power Skinning Algorithm.

Figure 14 :

 14 Figure 14: Motion LoD.

d

  = length(charactereye) T = α (t) = d-dnear d f ar-dnear ⇐⇒ α (t) = te-t te-ts where t ∈ [ts,te] and α(t) ∈ [0, 1]

Figure 15 :

 15 Figure 15: Benefit of skeletal Lod for crowds.

Figure 16 :

 16 Figure 16: Benefit of mesh Lod for a character.

Table 1 :

 1 BindGraphMesh and decimation performances.

		triangles t initBGM t 10 * decim
	25	982	142 ms	1 ms
	4	169	4 ms	< 1 ms
	21	7186	7.725 s	9 ms
	24	832	105 ms	2 ms
	25	3928	2.230 s	4 ms
	21	7186	8.848 s	9 ms
	24	6656	11.6 s	12 ms

c The Eurographics Association 2008.