David Prévost
email: david.prevost@lurpa.ens-cachan.fr

Sylvain Lavernhe
email: sylvain.lavernhe@lurpa.ens-cachan.fr

Claire Lartigue
email: claire.lartigue@lurpa.ens-cachan.fr

Didier Dumur
email: didier.dumur@supelec.fr

Feed drive modelling for the simulation of tool path tracking in multi-axis High Speed Machining

Keywords: Simulation, High speed machining, Multi axis machining, Follow-up, Identification, Generalized Predictive Control

Within the context of High Speed Machining, it is essential to manage the trajectory generation to achieve both high surface quality and high productivity. As feed drives are one part of the set Machine tool -Numerical Controller, it is necessary to improve their performances to optimize feed drive dynamics during trajectory follow up. Hence, this paper deals with the modelling of the feed drive in the case of multi axis machining. This model can be used for the simulation of axis dynamics and tool-path tracking to tune parameters and optimize new frameworks of command strategies. A procedure of identification based on modern NC capabilities is presented and applied to industrial HSM centres. Efficiency of this modelling is assessed by experimental verifications on various representative trajectories. After implementing a Generalized Predictive Control, reliable simulations are performed thanks to the model. These simulations can then be used to tune parameters of this new framework according to the tool-path geometry.

Introduction

High-Speed Machining (HSM) is now a standard process for the machining of sculptured surfaces as it allows high productivity combined with a good geometrical and dimensional accuracy [START_REF] Krajnik | Modern machining of die and mold tools[END_REF]. However, the numerous stages of the process, from the tool trajectory computation to its execution on the machine tool, involves a large number of parameters, data transformations and information exchanges, which may affect global productivity and final part quality [START_REF] Schmitz | Case study: A comparison of error sources in high-speed milling[END_REF]. Generally, the process is divided in three main stages (Figure 1): the CAM stage, the NC trajectory generation stage, and the cutting process itself [START_REF] Tournier | Iso-scallop tool path generation in 5-axis milling[END_REF].

Figure 1 Process stages

The CAM stage transforms the CAD model into a trajectory according to the machining strategy giving the CL-file. Sources of errors are mainly due to the calculation methods of tool positioning, and to the approximation of the CAD model by a set of points. The use of polynomial trajectory limits such errors (Lartigue, 2004).

Starting from the tool trajectory defined in the CL-file, the second stage carries out axis displacements. This stage calculates of the axis command followed by the trajectory tracking by the feed drives. During the last stage, the actual tool movement relatively to the part surface is executed on the machine tool; the cutting process thus generates the final machined surfaces. Therefore, the second stage is a key stage as it ensures the link between the numerical domain of the tool trajectory calculation and the physical domain of tool path execution by the axes. As such, it deserves a detailed description.

Dugas proposes a description of the process through three main levels defined in Figure 2 [START_REF] Dugas | Simulation d'usinages de formes complexes[END_REF]. The level 1 (or numerical level) interprets in real time the program lines of the CL-file to generate axis commands. In this level, the NC unit calculates the position and velocity set points based on the kinematics limits of each axis. The level 2 (or Feed drive level) corresponds to axis cards and drives which perform axis commands, after the conversion of analogue data into digital data. The level 3 (or kinematics level) consists of the mechanical structure of the machine tool. The mechanical chain creates the relative movement between the tool and the part from the axis commands previously calculated. This last level is the entry point of the cutting stage.

Actual joint trajectories

… …

…

Within the context of virtual machine tool, the prediction of the global process behaviour relies on the modelling of each level [START_REF] Altintas | Virtual Machine Tool[END_REF]. Most authors focus on a specific part of the modelling. In a previous work, a first model of trajectory execution has been proposed (Lavernhe et al., 2008a[START_REF] Lavernhe | Optimization of 5-axis high-speed machining using a surface based approach[END_REF]. By integrating kinematics performance, this model allows the simulation of the actual tool path generation during machining by determination of position, velocity, acceleration and jerk of each axis. Results have in particular highlighted that the step of interpolation (level 1) is source of high decreases of the velocity along the tool path, and geometrical approximations involving deviations on the machined part. However, in order to improve the simulation with respect to the real behaviour, a modelling of the feed drives (level 2) separated from the modelling of the interpolator is necessary. It would thus be possible to predict from the set points calculated at the first level the actual tracking of the trajectory resulting from the feed drives and, consequently, the analysis of the tracking and contour errors [START_REF] Ramesh | Tracking and contour error control in CNC servo systems[END_REF].

The paper deals with the modelling of the feed drives (level 2) with the aim of predicting actual axis displacements. With such a model, it becomes possible to predict the final geometrical error as the sum of the errors related to each level.

Many studies for modelling of feed drives have been carried out in the literature [START_REF] Whalley | Hybrid modelling of machine tool axis drives[END_REF], Yeung et al. 2006b[START_REF] Barre | Dynamic behaviour modelling and control of agile machine tool[END_REF]. Most of these works are performed on open CNC architectures or on specific machine tools (robots or agile machines), and in most cases considering one or two translational axes only. In the paper, the axis drive model proposed is deliberately simple (but realistic) to be applied to any type of axis (translation or rotation) and whatever the industrial high-speed machine tool. The model is developed for the context of HSM, considering finishing conditions so that disturbances due to cutting loads can be neglected. The model is assessed through various examples involving 2D or 3D trajectories. The result is an efficient model of the complete servo system, including existing axis controllers on industrial machining centres. Moreover, as the model proposed is generic, the implementation of various command strategies such as Generalized Predictive Control (GPC) is possible. To illustrate this aspect, the design and tuning of a Generalized Predictive Control (GPC) structure is proposed and implemented for the position loop of each axis drive.

Dynamic model

Common feedback structures used in industrial CNC drives are often based on a cascaded control structure. A typical cascaded structure (Figure 3) consists of several nested loops with dynamics becoming faster from the external to the inner loops.

Numeric axis control

The current inner loop aims at controlling the motor torque as the torque is directly proportional to the armature current. Within the cascaded framework, its high bandwidth enables fast disturbances rejection. The current loop is closed using a proportional integral (PI) control. The velocity central loop is also closed using a PI controller. The position loop is the most external one with a larger sampling period. This last loop is closed using a proportional gain. Velocity and torque feedforward actions are added to the structure with the aim of reducing the tracking error resulting from the lack of integrative action in the position controller. The structure defined above is commonly used in industrial CNC drives with the major drawback that the structure is not open, essentially for historical and robustness reasons.

The proposed model relies on the feed drive structure presented in Figure 3. Axes are modelled in an independent way, neglecting the existing coupling interactions. Some studies focus on these phenomena but they concern specific machines such as welding robots [START_REF] Ramesh | Tracking and contour error control in CNC servo systems[END_REF] or agile structure machines [START_REF] Barre | Dynamic behaviour modelling and control of agile machine tool[END_REF]. Besides, the model assumes the joints to be perfect and the components are supposed to be non deformable rigid bodies. The complete and generic model is built for translational axes as well as rotational ones. The implementation of the model is performed in the Matlab/Simulink™ environment (see Appendix 1). Depending on the machine architecture, some axes can be affected by gravity. This effect is taken into account by adding a resistant load (force or torque) to the motor. This resistant load is constant for a vertical axis whereas it depends on the axis angular position for a tilting table. Moreover, system of gravity compensation may exist. It must be modelled depending on technological solutions chosen by the manufacturer. Considering a context of finishing operation, disturbances due to cutting forces are neglected compared to dynamical effects (inertias, masses).

As the model must be generic and simple, some physical linked to the mechanical behaviour such as stiffness of the joints, backlashes, vibrations, are not modelled at this step (Siemens®, 2002;[START_REF] Whalley | Hybrid modelling of machine tool axis drives[END_REF]Yeung et al., 2006a). Numerical filters and numerical compensations are aggregated into three adjustment delay parameters, (α, β and γ) set in both the position and the feedforward loops. For each axis, the input of the model is the setpoint calculated by the interpolator and transmitted to the axis cards, so-called Position Setpoint at Entrance of the Controller (PSEC). The choice of this setpoint is adopted in order to be free from numerical operations that are specific to each manufacturer [START_REF] Bloch | Innovative feed rate optimisation technique[END_REF]Siemens®, 2002). The main output of the model is the simulated position (SP). Other inner variables of the model, such as the Simulated Velocity (SV) or the Simulated Motor Current (SMC) can be collected during simulation. In the next section, the identification step of the model parameters is detailed.

This model could be declined for different types of machine tool, with different motorization technologies. In the paper, the identification method is applied to one machine tool, which is a 5 axis Mikron UCP 710 HSM centre located in the LURPA lab, equipped with a RRTTT structure, rotary synchronous motors and Siemens Sinumerik 840D power line. Only the application to this 5-axis milling centre with a classical motorization is developed in the following sections.

Model parameter identification

Identification process

Generally, all the feedback controller parameters, such as the PI controllers, the proportional gain or the sampling periods are given by the CNC documentation. The parameters to be identified are thus the equivalent inertia, the friction coefficients and the feedforward compensations. In literature, many techniques for parameter identification exist. Some of them rely on frequency analysis by disconnecting the servo loops [START_REF] Erkorkmaz | High speed CNC system design. Part II: modelling and identification of feed drives[END_REF]Yeung et al., 2006a). Global techniques are also proposed so that the overall closed-loop dynamics is determined by running a standardized G-code and capturing the commanded and measured positions [START_REF] Erkorkmaz | Rapid identification technique for virtual CNC drives[END_REF]. However, these techniques can only be easily performed for open architecture machines, which is not the case for industrial machining centres. Modern CNC systems offer the possibility to capture various variables in real-time during the trajectory execution, such as current, position, velocity and acceleration. This functionality is used during the identification step but also for the determination of the feedforward actions (e.g. Torque Feed Forward Set point TFFWS and Velocity Feed Forward Set point, VFFWS).

Motor modelling

Motors, either linear or rotational, can be modelled by direct current motors, which are sufficiently precise as suggested in [START_REF] Boldea | Linear electric actuators and generators[END_REF]. Typical equations for direct current motors are given in (1). They are transformed in the Laplace domain to be implemented in the proposed model (see Appendix 1).

( ) . ( ) ( ( )) ( ) . ( ) . ( )
( ) . ( ) ( ( )) . ( ) ( )                 em mt m eq m r E t K t d i t U t R i t L E t dt C t K i t dt J C t C t dt (1)
E is the counter electromotive force, Ω m the instantaneous rotation velocity, K e the constant of counter electromotive force. U is the voltage across armature terminals, R is the electric resistance, L is the inductance, i the current. C m stands for the motor torque, and K t is the torque constant. J eq is the equivalent inertia projected on the motor axis and C r gathers all resistant torques.

Friction modelling

In axis control, in the absence of cutting forces, the dominant source of disturbance is friction. Several studies lead to the modelling of friction through relationships between friction, position, speed and temperature of rigid bodies [START_REF] Kikuuwe | Fixed-step friction simulation: from classical Coulomb model to modern continuous models[END_REF][START_REF] Olsson | Friction models and friction compensation[END_REF]. In this paper, a friction model depending only on axis velocity is chosen, as suggested by [START_REF] Olsson | Friction models and friction compensation[END_REF] and (Mennon et al., 1999). The friction model is the combination of both viscous and Coulomb friction. When acceleration is zero, the motor torque is equal to the resistant torque. Considering equations (1), this yields to:

( ) ( ) . ( )  m r t C t C t K i t (2)
If the resistant torque is only caused by friction, the equivalent current called i friction can be identified through NC recordings. Indeed, it is possible to evaluate the actual friction law by recording the current i friction for various displacements with constant velocities. Trials performed on the machine tool show that i friction is independent of the axis position and only depends on the velocity (Figure 4).

                   b V d V friction b V d V friction friction i V a e c e if V i V a e c e if V i i i if V
(3)

Equivalent inertia and masses

Concerning the motor parameters, electric parameters can be found in manufacturer's documentation. An experimental identification is performed for mechanical parameters, which account for embedded masses and other inertias.

Equivalent inertia relative to the motor shaft is calculated from a series of current measurements with zero acceleration displacements according to equation ( 5). J eq is then calculated for each time step and an average is computed on each time interval for which acceleration is non zero.

.

   tr eq m K i C J d dt (5)
As previously suggested, gravity is taken into account, especially for vertical axes and tilting axes, by adding resistant torques. These ones are easily determined with a measurement of the current on static positions, when there is neither velocity nor acceleration.

Feedforward and shift parameters

Feedforward is a common technique used to cancel dragging differences on the position [START_REF] Lambrechts | Trajectory planning and feedforward design for electromechanical motion systems[END_REF]. Two types of feedforward are commonly implemented and will be modelled bellow: the torque feedforward and the velocity feedforward, which can be disabled or used separately or together. Therefore, constants associated to the feedforward compensations can be evaluated from the measurement of the feedforward setpoints TFFWS and VFFWS. The ratios between the first and second derivatives of the Position Setpoint at Entrance of the Controller (PSEC) and VFFWS or TFFWS respectively are supposed to be constant:

( ( )) ( ) .  d P SE C t V F F W S t V F F W dt (6) 2 2 ( ( )) ( ) .  d P SE C t T F F W S t T F F W dt (7)
Derivatives are computed according to the Euler approximation (with z -1 the backward shift operator and Tps the sampling rate):

1 2 1 2 22 ( ( )) 1 ( ) ( ) . ( ) ( ( )) 1 2. ( ) ( ) . ( )           d x t z y t Y z X z dt T ps d x t z z u t U z X z dt T ps (8)
Then, a linear regression is applied to measurements to verify that feedforward terms are constant. Following the identification stage, an adjustment of the parameters is performed. To do this, three parameters (α, β, γ) are introduced (see Appendix 1), which stand for a delay in the feed drive structure. They are optimised using a leastsquare minimisation of the deviations between simulated and measured positions. To begin with, the first parameter α located in the position loop is identified without feedforward actions. Next, β is identified in the same way with only the velocity feedforward and finally, the last one γ when all anticipations are activated. Values of adjustment parameters are given in Appendix 2.

This section details the model assessment for the 5-axis milling centre Mikron UCP710 performed thanks to trials. Different types of trajectories are tested, involving one or several axes. Data are directly recorded via the NC oscilloscope at the position loop sampling period. The measured signals are the position setpoints (PSEC), the actual position recorded through the encoders (mm), the current (A) and the velocity (m/min).Collected data are directly confronted to values issued from simulations.

One axis validation

Testing one-axis permits to check independently each axis while being free from coupling problems. Trials are performed for translational axes as well as rotational ones. Whatever the axis, all the experiments lead to similar results. Therefore, only tests concerning the X translational axis and the C rotational axis are presented in the paper. Considering a segment defined from two position values, several tests are conducted:

 Case 1: the programmed velocity is constant during the displacement (Figure 6)  Case 2: the programmed velocity is V1 for a segment portion, and is V2 for the other part, with V1>V2 (Figure 7)

 Case 3: the segment is travelled in back and forth at the same programmed velocity.

The case 1 is illustrated for the X axis in Figure 6. Simulated values are closed to the measured ones as deviations do not exceed 4 µm except during critical phases of acceleration and deceleration. In such cases, differences can reach 6 µm. Simulated and measured values for the velocity and current are almost identical. Hence, values are small enough to assess the model of each translational axis independently, with or without feedforward actions.

Case 2 is illustrated for the rotational C axis in Figure 7. Practically, this test consists in covering the angular interval [0° , 130°] at the speed of 18 rpm and the interval [130° , 210°] at 6 rpm. Discrepancies collected between angular positions measured and simulated are smaller than 6.10 -3 degrees for this axis, which is rather a good result. Once more, the maximum deviations are obtained during deceleration and acceleration phases. Hence, whatever the axis nature (translational or rotational), and the motorization (synchronous rotational or synchronous linear motors), differences between measurements and simulations are sufficiently small to assess the model for each axis independently solicited. The accuracy of the simulation (globally less than 0.005 mm or 0.005°) seems to be enough to predict tracking deviations.

Multi axis tests

These tests aim at the model assessment for trajectories involving several axes at the same time. Moreover, trajectories present some geometrical singularities implying significant tracking deviations in particular for high velocities.

The four trajectories tested are presented in Figure 8. For each case, differences between simulated and measured values are estimated (right part of the figures). Deviations between simulations and measurements never exceed 0.01 mm for translational axes, respectively 0.01° for rotational axes; deviations are below 5 µm (resp. 0.005°) for all the trials. However, the deviations reach 12 µm for the 3D curve. Such differences are due to chatters that are discernible during acceleration or deceleration phases when the trajectory is executed on the machine tool.

Differences observed between simulated and measured data are likely due to the difficulty of modelling the specific numerical and software NC treatments. Indeed, such elements are considered through the adjustment parameters (a, b, g) which are not so simple to precisely identify. Furthermore, accounting for axis coupling or backlash compensations could also improve the model precision.

Nevertheless, the comparisons proposed in Figure 8 enhance the good behaviour of the model to simulate trajectory tracking whatever the trajectory executed. Moreover, the model allows the prediction of tracking deviations that can impact the geometrical quality of the parts [START_REF] Prévost | Feed drive simulation for the prediction of the tool path follow up in High Speed Machining[END_REF].

Figure 8 Multi axis tests

In the following, one application using the model presented above is studied. It deals with the implementation and the tuning of advanced controllers, and more particularly predictive controllers based on a GPC structure.

Towards an optimized RST Generalized Predictive Control

This section presents the implementation of an advanced controller structure based on the GPC algorithm. Among various strategies of advanced controllers, previous research works in the field of machine tool in particular show that GPC can provide better tracking performance compared with classical cascaded structures (** [START_REF] Susanu | Improving performance of machine tools with predictive axis controllers within an open architecture Framework[END_REF].

However, these works also highlight that the choice of the required tuning parameters is a crucial point to obtain good tool path tracking performance, which may depends not only on the behavior of the controlled loops, but also on the kind of trajectory issued from the interpolation module. On the other hand, implementing predictive strategies is performed with the help of a prediction model of the system, which is usually selected as simple as possible. These reasons are strong motivations in favor of the model developed in this paper, which will be very helpful in particular from a tuning point of view. This aspect will be further examined in the specific case of the Mikron UCP 710 machining centre.

A first paragraph summarizes the basic ideas of predictive control, leading to the elaboration of the predictive controller under the RST structure. A second paragraph then details the implementation of such a strategy in the case of trajectories generated through circular interpolation, emphasizing the need for optimization of the parameters depending on the type of trajectory.

GPC strategy and implementation of the new structure

Based on Model Predictive Control philosophy, the Generalized Predictive Control strategy relies on four important ideas reproducing the basic decisional of the human behaviour [START_REF] Rossiter | Model -Based Predictive Control[END_REF], creation of an anticipative effect by exploiting the trajectory to be followed in the future, definition of a numerical representation of the process as the prediction model, minimization over a finite horizon of an objective function elaborating the best control action, timedomain displacement of predictions known as receding horizon strategy. These basic steps can be translated in equations and algorithms which are summarized in [START_REF] Susanu | Hierarchical Predictive Control within an Open Architecture Virtual Machine Tool[END_REF], leading to a controller finally given in the socalled RST formalism. These polynomials R, S and T are completely determined by algorithms that will not be detailed here and whose coefficients derived from the resolution of Diophantine equations. As a result, the tuning of four parameters N 1 , N 2 , N u and λ (see Table 1) is required prior to the elaboration of the three polynomials R, S and T, which are then computed offline and uniquely defined. With this structure, the T polynomial has a non causal form, thus ensuring an anticipative effect in closed-loop. Consequently, there is no need anymore for implementing feedforward actions which are difficult to tune, as usually done in classical NC. Finally the RST controller is implemented as shown in Figure 9. In the following test, the RST formalism is implemented in a cascaded version, only for the position loop of each axis drive.

Figure 9

The RST polynomial structure Feedback is thus carried out under a numerical RST structure, only requiring the computation of a simple difference equation:

) ( ) ( ) ( ) ( ) ( ) ( * 1 1 * 1 t x q T t x q R t q S        (9)
where x* is the input, x the output and * the velocity reference signal.

Remark 1: Using this RST formalism to structure the axis control law represents an important advantage since this form appears to be very generic. Indeed, a large number of numerical controllers can be structured in the same way, from basic PID to more advanced ones. Within the framework of open architecture, axis controllers can thus be seen as plug and play modules, ensuring interchangeability features.

Remark 2: Another interest of this generic structure is the resulting reduced computation time as the controller polynomials are in most cases of small degrees. This aspect allows the implementation of this strategy even with short sampling times, as for instance in the case of very high speed machining.

Based on the feed drive model developed in the previous sections, the following section presents the results of simulation with the predictive RST control structure on the position loop, considering a tuning of the four parameters N 1 , N 2 , N u and λ ensuring good stability and robustness properties (** [START_REF] Susanu | Improving performance of machine tools with predictive axis controllers within an open architecture Framework[END_REF]. Comparisons are made between classical control architecture with feedforward terms and predictive RST control structure for two types of trajectories.

RST simulations and comparisons with the classical control architecture

The first case consists in plane circular trajectories, so-called ballbar test, largely used to analyse and measure tracking errors and geometrical errors of the machine tool. A circular trajectory of radius 150mm is first considered, with a programmed feedrate of 15m/min. The trajectory is interpolated using circular interpolation (G2/G3) and with velocity and torque feedforward actions on. Figure 10 represents tracking errors with classical cascaded control structure and feedforward actions (left) and predictive RST control structure (right).

Figure 10 Tracking error comparison between cascaded and RST structurecircular trajectory

It can be clearly seen that the 50mm tracking error observed with the cascaded structure becomes almost zero with the RST structure. As a first conclusion, the chosen parameters tuning of the GPC structure greatly improves the axes behaviours. The fact that the industrial NC module is completely closed does not enable to validate the predictive law through real experiments. In that sense, the elaboration of a realistic and simple model as proposed in the paper is a major advantage which helps validating new advanced controllers.

The second case, for which trajectories present high tangency discontinuities, is now studied. For example, two trajectories in the X-Y plane represented on Figure 11 are tested. The geometries chosen here are sharp corners (45° and 90°), which generally induce high deviations on machined parts with high programmed feedrates. The programmed feedrate speed for each case is 10m/min. Figure 11 presents the behaviour of the simulated trajectories, obtained for the predictive controller with the same tuning as for the first case. A decrease of the gap between simulated position with the new control strategy (GPC) and the position setpoint is illustrated on the right handside of the figure. Unfortunately, this is not the case for the left handside of this figure. This case highlights the need for optimization of the tuning of the four parameters N 1 , N 2 , N u and λ. This can be realized using the model developed in this paper. Indeed, the virtual environment developed around the feed drive model is a useful tool to select through simulation the adequate tuning for the set of parameters, in particular in function of the kind of trajectory that has to be followed. One perspective is thus to find an optimisation technique and propose such a technique in function of the tool path typology.

Conclusion

This paper focuses on a model of the feed drives of machine tool. The proposed model is deliberately simple to be applied to milling centers equipped with modern industrial CNC to simulate tracking errors. The model accounts for the main influent components on the dynamical axis behaviour: inertia, friction, and feedforward compensations. Complex phenomena as for specific NC treatments are taken into account through adjustment parameters. An original procedure of parameter identification is carried out thanks to the capture on fly of various variables using the oscilloscope of the NC unit. Therefore, the virtual environment developed around the feed drive model gives a useful tool to simulate tracking errors, but also to investigate various control structures.

The model is implemented with success for a 5-axis milling center. Deviations observed between the measured displacements and the simulated ones relatively low assess the relevancy of the model. To enhance the benefit of such a model various types of trajectories are tested. Results clearly emphasize the good behaviour of the model to simulate trajectory tracking whatever the trajectory executed, even for 3D trajectories. The model turns out to be useful for the prediction of tracking deviations that can impact the geometrical quality of the parts.

A second major interest of the model is its adaptability to various control structures. In this direction a Generalized Predictive Control is implemented into the model to design a new framework of control. From the control point of view, the elaboration of a virtual environment based on a simple but realistic model of the feed drives presents at least two important advantages. First, it allows testing and validating through realistic simulations new advanced control architectures, which would not have been possible up to now on industrial machining centers due to the closed architecture of the NC modules. Then, it allows through several simulations selecting the best set of tuning parameters to provide stability and robustness in relation with a specific trajectory.

Figure 2

 2 Figure 2 Global architecture of a NC machine tool

Figure 4

 4 Figure 4 Current and speed measurement (X axis)

Figure 5

 5 Figure 5 Friction law and experimental results (X axis)

Figure 6 Figure 7

 67 Figure 6 X axis translation test

  Figure 8.a and 8.b report results concerning trajectories executed in the X-Y plane, with a programmed feedrate Vf =10 m/min. Figure 8.c stands for tool orientation change in the part coordinate system considering a trajectory defined in the A-C plane. Finally, in Figure 8.d, the trajectory is a 3D B-Spline curve.

Figure 11

 11 Figure 11 Trajectory comparison between cascaded and RST structuresharp corners

Figure 3 Axis control structure

											Numeric axis control
								Torque Torque
							Feed Forward Feed Forward
					Velocity Velocity
				Feed Forward Feed Forward
	Interpolator Interpolator	x* y* y* x*	+ +	--	ε p ε p	Position Controller Position Controller	+ + + +	--	ε v ε v	Velocity Controller Velocity Controller	+ + + +	--	ε a ε a	Current Controller Current Controller	u u	Axis motor and load Axis motor and load	i i	Ω x x Ω
						Axis control Axis control				Axis motor and load Axis motor and load		y y
		z* z*
						Axis control Axis control				Axis motor and load Axis motor and load		z z

Table 1

 1 RST tuning parameters

	N 1	Minimum output prediction horizon
	N 2	Maximum output prediction horizon
	N u	Control horizon
	λ	Weighting factor

Acknowledgements

Authors would like to thank AIP Primeca Ile de France who have graciously allowed to use of the 3x DMG HSC75 milling center. This work was carried out within the Manufacturing 21 working group, which comprises 18 French laboratories. The topics approached are:

-modelling, -virtual machining, -emergence of new manufacturing methods.


Zhengdong Zhou

Pascal Haigron

Huazhong Shu

Wenxue Yu

Oscar Acosta

Jean-Pierre Manens

Antoine Lucas

Liming Luo

A Study of Intravascular Brachytherapy Treatment Planning in Peripheral Arteries

A two-step method to determine the seed parameters for the planning of peripheral intravascular brachytherapy, which took the actual vessel geometry into consideration, was developed. Firstly, the characteristics of the vessel geometry was obtained by using active navigation; Then a method combined genetic algorithm with BFGS algorithm was applied to optimize the number of seeds, and the parameters associated with each seed such as the position, and the dwell time. Application of the method to a phantom model and three animal models of stenosis shows that promising result could be obtained, and the planning of peripheral intravascular brachytherapy should take actual vessel geometry into consideration.

I. INTRODUCTION

Intravascular Brachytherapy (IVB) constitutes a new therapeutic solution to avoid iterative redilations practiced after minimally invasive treatment of peripheral arterial stenoses. Following percutaneous transluminal angioplasty (PTA) procedure, IVB aims at preventing intimal hyperplasia restenosis in a context which does not tolerate iatrogenic complications.

Although its principle is known, the application of IVB to peripheral arteries by means of a high dose rate (HDR) remote afterloader ( 192 Ir seed) remains at a preliminary stage. Peripheral arteries present some specific features compared to coronary arteries as for their shape, geometrical characteristics and the presence of surrounding anatomical structures. Conventional treatment planning still assumes a straight line vessel and a constant diameter circular crosssection. In order to use IVB on peripheral arteries in clinical practice, several issues have to be considered: (i) definition of the protocol, (ii) anatomical feature based dosimetric planning, (iii) implementation of the treatment plan, and (iv) evaluation.

This work deals with the optimization of the treatment plan for IVB in peripheral arteries. The objective is both to evaluate the difference between the straight line hypotheses and actual vessel geometry, and to contribute to the definition and the specification to a new optimized procedure. In the proposed approach, the planning process is divided into two steps: (i) characterization of the vessel geometrical parameters by active navigation, and (ii) optimization of the seed parameters using simulated annealing and BFGS algorithm. The first step was presented in a previous work [START_REF] Haigron | Planning of intravascular brachytherapy based on virtual exploratory navigation in X-ray CT images[END_REF].

We focus here on the optimization process in which the number of seeds, the position and dwell time associated with each seed are taken into account.

II. RELATED WORKS

The results of researches on IVB showed that less renarrowing of arteries and a reduced need for further treatment at the same site in patients with in-stent restenosis who received radiation along with angioplasty. And now intravascular brachytherapy has become a viable treatment option in preventing restenosis.

Nevertheless, in most current planning procedures of intravascular brachytherapy, the lesion path is supposed to be straight line and the dose distribution is specified with respect to the seed and the mean diameter of the artery. In order to make intravascular brachytherapy procedures more secure and reliable, recent works, on intracoronary brachytherapy notably, have shown that the morphological characteristics of the vessel should be taken into account in the planning process, because the doses closely depend on the geometry of vascular [START_REF] Nath | Intravascular brachytherapy physics: Report of the AAPM Radiation Therapy Committee Task Group N 60[END_REF][START_REF] Wahle | Effects of vessel geometry and catheter position on dose delivery in intracoronary brachytherapy[END_REF]. One solution adopted in the case of coronary arteries is using IVUS to detect the inner surface of the vascular structure [START_REF] Nath | Intravascular brachytherapy physics: Report of the AAPM Radiation Therapy Committee Task Group N 60[END_REF][START_REF] Roelandt | Three-dimensional reconstruction of intracoronary ultrasound images: rationale, approaches, problems and directions[END_REF][START_REF] Wahle | Estimating the actual dose delivered by intravascular coronary brachytherapy using geometrically correct 3-D modelling[END_REF]. In the case of peripheral arteries, other solutions such as CT imaging can be used to obtain a better accuracy in 3-D boundary detection and description. We investigated the use of a solution based on active navigation, a virtual angioscopy like process which does not require a pre-processing of the patient-specific CT data. The general framework of active navigation, whose basic functions have already been reported in the literature [START_REF] Haigron | 11Depth-map-based scene analysis for active navigation in virtual angioscopy[END_REF], is summarized in the following.

Virtual image computation is performed according to a ray casting scheme integrating analysis functions. From the depth map associated to the virtual image, the scene analysis process is based on the detection of the depth maxima charaterizing the vascular branches and border points featuring, in a bifurcation area, the frontier between two different regions associated to two single branches respectively. The depth maxima are labelled and matched along the succesive positions of the virtual endoscope, so that, given a starting point and a target point, a particular branch can be followed. The viewing direction of the virtual sensor is oriented towards the followed depth maximum whereas the displacement orientation is fixed initially by the center of a set of surface points selected in a single branch part.

The geometrical model of the virtual angioscope has been extended to include a quantification plane which allow it to construct, during the navigation, a precise surface description of the vascular structure [START_REF] Acosta | Geometrical description of anatomical structures based on virtual exploratory navigation[END_REF].

III. SEED MODEL

To determine the parameters that may influence the treatment plan, it is first necessary to choose an appropriate formulae for computing the dose distribution at an arbitrary point. The standard Nucletron 192 Ir cylindrical seed we considered is 3.5 mm long and 0.6 mm in diameter, encapsulated in a 5.0 mm long × 1.1 mm diameter stainless steel housing.

The computation of dose distribution often uses the basic dosimetry data, including dose rate constants, radial dose functions, and anisotropy functions. According to the formalism proposed and recommended by the AAPM in TG43 and TG60 [START_REF] Nath | Intravascular brachytherapy physics: Report of the AAPM Radiation Therapy Committee Task Group N 60[END_REF][START_REF] Nath | Dosimetry of intersticial brachytherapy sources: Recommandations of the AAPM Radiation Therapy Committee Task Group N 43[END_REF], the dose distribution, for a cylindrically symmetric source, can be described in terms of a polar coordinate system with its origin at the source center where r is the distance to the point of interest and θ is the angle with respect to the long axis of the source. The dose rate, D(r, θ), at point (r, θ), is defined as follows:

D(r, θ) = S k [G(r, θ)/G(r 0 , θ 0 )]g(r)F (r, θ) (1)
where S k is the air kerma strength of the source, Λ is the dose rate constant, G(r, θ) is the geometry factor, g(r) is the radial dose function, and F (r, θ) is the anisotropy function. The reference point (r 0 , θ 0 ) is chosen to lie on the transverse bisector of the source at a distance of 1cm from its center, i.e., r 0 = 1 cm and θ 0 = π/2.

When the number of sources used in a treatment is N , the dose rate received at a point P is the weighting sum of D(r, θ):

I(P ) = N j=1 w j D(r j , θ j ) (2)
where w j corresponds to the irradiation time of the jth source.

IV. IVB PLANNING OPTIMIZATION

In order to fit the calculated dose distribution to the dose prescribed according to the 3D surface description of the vessel wall, and to minimize the damage to critical or normal tissues, a quadratic, least-square objective function was used to optimize the irradiation parameters. Let X denotes the set of parameters to be optimized, we have:

minf (X ) = α NT i=1 (D T -I(p i )) 2 + β NN i=1 δ(I(P i ) -D N )(I(P i ) -D N ) 2 (3) with δ = 1 if x > 0 0 otherwise subject to w j ≥ 0
where N T is the number of 3D surface sampled points in the lesion area; N N is the number of 3D surface sampled points distributed on the normal tissues; D T is the prescribed dose at the point P i of the lesion surface; D N is the prescribed dose at the point P i of the normal tissues; I(P i ) is the calculated dose at point P i . The coefficients α and β are the degree of desirability weighting factors for target points and normal or critical points, respectively.

Many algorithms are available to solve the optimization problem (3), most of them are stochastic methods. Here we firstly replaced the parameter w j with y 2 to eliminate the problem of solutions with negative dwell time, thus the problem (3) became an unconstrained optimization problem; then a method combined simulated annealing algorithm with BFGS algorithm was used to optimize the number of seeds, the dwell position as well as the dwell time associated with each seed.

In this paper, BFGS algorithm was used to determine the dwell times for a given number of seeds and corresponding dwell positions, while simulated annealing algorithm was used to determine the positions of the seeds, the parameters were set as follows; the initial temperature was set to 1000, the cooling factor was set to 0.85, the final temperature was set to 1.0, and the number of iterations per fixed temperature was set to 100.

The range of the number of seeds depends on both the length of plan path, given here by the vessel centreline, and the seed length. Let L s be the seed length, and L p the plan path length, the number of seeds is assumed to vary from L p /(1.2 × L s ) to L p /(0.8 × L s ). For each given number of seeds, Eq. ( 3) was solved by the combined method to determine the corresponding optimal solution. Among them, the one fits best to the prescribed dose distribution was selected as the final result.

V. EXPERIMENTAL RESULTS

In order to evaluate our methodology, we considered four data sets acquired from helical CT scanner, one from a phantom model, one from pre-PTA iliac artery and two from post-PTA sheep iliac artery. Several kinds of parameters intervene in the proposed approach. The aim is to evaluate the difference of the result between the straight line assumption and the actual vessel geometry. The criteria considered to quantify the results of the planning process is the dose surface histogram(DSH) which represents the dose homogeneity over the treatment area. The geometrical description of vascular inner surface was detected by active navigation. The activity of the seed was set to the maximum possible value (10ci or 370Gb).

A dose of 14Gy has been specified at a distance of 2 mm to each considered surface points (radius + 2 mm). Preliminary results are reported in figure 1(a-l), (a), (b), (c), (d) are the automatic reconstruction of the vessel inner surfaces and their center lines for phantom model, pre-PTA iliac artery and two from post-PTA sheep iliac artery respectively; (e), (g), (i), (k) are the optimized DSHs corresponding to the above four data set, taking actual vessel geometry into consideration; (f), (h), (j), (l) are the optimized DSHs corresponding to the above four data set, with straight line assumption. The results list in figure 1 show that DSHs taking actual vessel geometry into consideration are much better than those with straight line assumption.

VI. CONCLUSION AND FUTURE WORKS

This study demonstrates that the method can give promising result, and the planning of peripheral intravascular brachytherapy should take actual vessel geometry into consideration. Future work will take into account the currently available afterloader parameters.

Fig. 1 .

 1 Fig. 1. (a), (b), (c), (d) Automatic reconstruction of the inner surfaces of vessels and their center lines; (e), (g), (i), (k) optimized DSHs taking actual vessel geometry into consideration; (f), (h), (j), (l) optimized DSHs with straight line assumption.


Andrea Rueda

Oscar Acosta

Michel Couprie

Pierrick Bourgeat

Jurgen Fripp

Nicholas Dowson

Eduardo Romero

Olivier Salvado

Topology-corrected segmentation and local intensity estimates for improved partial volume classification of brain

Keywords: Brain tissue segmentation, Partial volume classification, Magnetic resonance imaging, Topology correction, Sulci detection, Cortical thickness estimation

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Accurate segmentation of Magnetic Resonance (MR) images into different brain tissues, namely gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF), can allow in-vivo quantification of structural modifications appearing during neurodegenerative diseases. However, MR-related artifacts, such as intensity inhomogeneity, noise and partial volume (PV) effects, can hamper the precision of this task. Inhomogeneities can be characterized by a low frequency multiplicative bias field and are mostly due to the sensitivity of the receiver coils and, in some cases, to non-homogeneous tissue MR properties. The noise is Rician distributed and it has be shown to strongly affect the tissue classification [START_REF] Van Leemput | A unifying framework for partial volume segmentation of brain MR images[END_REF]. Finally, PV effects appear when the size of anatomical features being imaged is comparable to the voxel size, causing blurring at the interfaces between tissues. In some cases, e.g. with opposed banks of GM in deep sulci, misclassification problems appears, affecting further processings such as cortical thickness estimation.

Topological operators and constraints have been widely used to correct and achieve accurate cortical tissue segmentations [START_REF] Ségonne | Active contours under topology control -genus preserving level sets[END_REF][START_REF] Bazin | Topology correction using fast marching methods and its application to brain segmentation[END_REF][START_REF] Han | Topology correction in brain cortex segmentation using a multiscale, graph-based algorithm[END_REF][START_REF] Kriegeskorte | An efficient algorithm for topologically correct segmentation of the cortical sheet in anatomical mr volumes[END_REF]. It has been assumed that the cerebral cortex is a folded sheet of GM built upon the WM, which would have the topology of a hollow sphere if the midline hemispheric connections were artificially removed. Due to MR artifacts, the segmentation process cannot guarantee this assumption, generating deviations from the true anatomy of the structures of interest. Proposed approaches that address this issue can be classified in two categories: methods that include topological constraints directly into the segmentation process, based on active contours [START_REF] Ségonne | Active contours under topology control -genus preserving level sets[END_REF], topology adaptive snakes [START_REF] Mcinerney | Topology Adaptive Deformable Surfaces for Medical Image Volume Segmentation[END_REF], digital topology models (Bazin andPham, 2005, 2007) or segmentation by registration to an atlas [START_REF] Kriegeskorte | An efficient algorithm for topologically correct segmentation of the cortical sheet in anatomical mr volumes[END_REF]; and retrospective techniques that correct the topology after the segmentation process [START_REF] Han | Topology correction in brain cortex segmentation using a multiscale, graph-based algorithm[END_REF]. Those approaches are focused on ameliorating the topology of the segmented tissues, working directly on a voxel or on a mesh (surface) space. Voxel-based methods operate directly on the volumetric tissue segmentations, by removing or adding voxels according to topological constraints. However, remotion or addition of a whole voxel in thin structures such as the GM may considerably modify the measure of thickness (ranging between ± 1 voxel) if any mechanism such as partial volume is not used to compensate for the structural modifications. In contrast, mesh-based techniques requires an initial 3D reconstruction (triangular mesh) of the volumetric segmentations. The approaches for segmentation and cortical thickness estimation operating directly with the surfaces, such as CLASP [START_REF] Kim | Automated 3D extraction and evaluation of the inner and outer cortical surfaces using a Laplacian map and partial volume effect classification[END_REF], BrainVISA [START_REF] Mangin | From 3D magnetic resonance images to structural representations of the cortex topography using topology preserving deformations[END_REF] or Freesurfer [START_REF] Dale | Cortical surface-based analysis I: Segmentation and surface reconstruction[END_REF][START_REF] Fischl | Cortical Surface-Based Analysis II: Inflation, Flattening, and a Surface-Based Coordinate System[END_REF][START_REF] Fischl | Measuring the thickness of the human cerebral cortex from magnetic resonance images[END_REF], incorporate mechanisms to pre-vent self-intersection of surfaces or topology correction, imposing also some smoothness constraints. Mesh-based approaches are however computationally more expensive, because of the needed additional reconstruction step.

Overall, after or during the mesh generation, most of the methods tackle the elimination of tunnels and handles [START_REF] Fischl | Cortical Surface-Based Analysis II: Inflation, Flattening, and a Surface-Based Coordinate System[END_REF][START_REF] Florent Ségonne | Geometrically accurate topologycorrection of cortical surfaces using nonseparating loops[END_REF][START_REF] Jaume | Open Topology: A Toolkit for Brain Isosurface Correction[END_REF][START_REF] Zhou | Topology repair of solid models using skeletons[END_REF].

On the other hand, PV estimation has received considerable attention in the last few years and different approaches have been proposed for classification and computation of fractional content [START_REF] Santago | Quantification of MR brain images by mixture density and partial volume modeling[END_REF][START_REF] Laidlaw | Partial-volume Bayesian classification of material mixtures in MR volume data using voxel histograms[END_REF][START_REF] Shattuck | Magnetic resonance image tissue classification using a partial volume model[END_REF][START_REF] Noe | Partial volume segmentation of cerebral MRI scans with mixture model clustering[END_REF][START_REF] Van Leemput | A unifying framework for partial volume segmentation of brain MR images[END_REF][START_REF] Tohka | Fast and robust parameter estimation for statistical partial volume models in brain MRI[END_REF][START_REF] Chiverton | Adaptive partial volume classification of mri data[END_REF]. Most techniques model voxel intensity as a linear combination of the intensity distributions of the possible tissue types within each voxel [START_REF] Choi | Partial volume tissue classification of multichannel magnetic resonance images -a mixel model[END_REF][START_REF] Noe | Partial volume segmentation of cerebral MRI scans with mixture model clustering[END_REF]. Computing the fractional content of voxels therefore requires both pure and mixed voxels to have been previously classified. [START_REF] Shattuck | Magnetic resonance image tissue classification using a partial volume model[END_REF] implemented a maximum a posteriori (MAP) classifier, which combined a tissue measurement model with a prior model of the local spatial interactions to obtain six tissue types: three pure and three mixed. The fractional content for the mixed voxels was calculated based on the global intensity mean of pure tissue types. [START_REF] Tohka | Fast and robust parameter estimation for statistical partial volume models in brain MRI[END_REF] proposed an algorithm which used statistical estimators, based on the MAP estimation [START_REF] Shattuck | Magnetic resonance image tissue classification using a partial volume model[END_REF]. Recently, [START_REF] Chiverton | Adaptive partial volume classification of mri data[END_REF] The approaches previously presented have been focused on solving either the PV estimation or the topology correction. Our contribution consists in demonstrating that better results and performance are obtained if both strategies are combined together with a spatial intensity variation modeling.

In this paper, we propose a new method aimed at improving both PV classification and fractional content computation, working at a voxel level in order to be accurate and computationally efficient. The improved classification is achieved by imposing topological constraints to the binary segmentation and thus detecting hidden mixed voxels in zones of tight sulci. The accurate fractional content estimation is attained by computing the fractional content as a linear relation between robust local intensity averages of pure tissue voxels. The spatially dependent averaging helps to overcome the problems of intensity inhomogeneity for a given tissue across the image.

In the next section we describe our methods, followed by experiments using simulated and real data. We also compare the results with other previously proposed methods. We demonstrated the utility of our approach by integrating the whole process to our voxel-based cortical thickness estimation pipeline.

Methods

The proposed strategy follows the steps depicted in Figure .1: Firstly, an initial classification of voxels into pure tissues WM, GM and CSF and mixed tissues WM/GM and GM/CSF is performed. Secondly, topology-constraints are introduced in the classification assuming that the GM is a continuous layer covering the WM. A topology preserving dilation of the WM over GM adds robustness to the delineation of mixed voxels GM/CSF in deep sulci.

Finally, the estimation of fractional content for mixed voxels is adaptively performed based on a local averaging of the pure tissue voxels.

INSERT FIGURE .1 HERE

Pure tissue segmentation

A first segmentation of pure brain tissues into GM, WM and CSF is performed based on an implementation of the expectation-maximisation (EM) segmentation method as in [START_REF] Van Leemput | Auto-mated model-based bias field correction of MR images of the brain[END_REF]. Here, the Colin atlas and associated priors are first affinely registered to the data using a robust block matching approach [START_REF] Ourselin | Reconstructing a 3D structure from serial histological sections[END_REF], followed by a diffeomorphic Demons non-rigid registration [START_REF] Vercauteren | Non-parametric diffeomorphic image resgistration with the demons algorithm[END_REF]. Probabilistic tissue maps associated with the atlas were used to initialize the EM segmentation and enforce spatial consistency throughout the segmentation. The probability density functions of the tissues are modelled with 6 Gaussians (WM, GM, CSF and 3 for non brain tissues, skull and background). Finally, hard segmentations are obtained after the EM segmentation by labelling each voxel with the most probable tissue.

Initial partial volume labelling

Using the hard segmentations, a first labelling of partial volume voxels are identified within the hard segmentations and along the interfaces of pure tissues. Three pure tissue classes and two mixture classes are considered Γ = {GM, CSF, WM, CSF/GM, GM/WM}. A maximum a posteriori classification (MAP) is made and labels the voxels as belonging to the set Γ. This procedure, relying on both intensity and spatial information, extends the method proposed by [START_REF] Shattuck | Magnetic resonance image tissue classification using a partial volume model[END_REF]), but we assume that each voxel contains at most two tissues [START_REF] Santago | Quantification of MR brain images by mixture density and partial volume modeling[END_REF], and PV classification is restricted to the region formed by a dilated GM region (radius 2) because only the cortical thickness is sought. To take into account dependency on the neighbouring tissue types, a Markov prior that models local spatial interactions was implemented using a Potts model in order to perform the labelling.

As in [START_REF] Shattuck | Magnetic resonance image tissue classification using a partial volume model[END_REF][START_REF] Tohka | Fast and robust parameter estimation for statistical partial volume models in brain MRI[END_REF][START_REF] Kim | Automated 3D extraction and evaluation of the inner and outer cortical surfaces using a Laplacian map and partial volume effect classification[END_REF], we use the Iterated Conditional Modes (ICM) algorithm as explained in [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF] to search for the optimal labelled image. According to this, every voxel is updated once per iteration until no label changes occur between iterations.

This model favors classification of contiguous regions of GM, WM and CSF and encourages configurations of voxels that make physical sense such as GM/CSF or GM/WM voxels adjacent to GM.

A topology preserving segmentation

After the MAP labelling, some of the sulci may be misdetected, as the intensity of buried PV GM/CSF voxels is close to that of the GM. In order to refine the segmentation and identify such buried GM/CSF voxels, we used a homotopic dilation of the consolidated WM = {WM, WM/GM} constrained by the GM, leading to a better delineation of deep sulci. To preserve this folds during dilation, the set WM is corrected first to assure that shares the topology of a filled sphere.

The homotopic transformations that we used are topology-preserving procedures that consist of sequentially deleting or adding single points (voxels) as described in [START_REF] Bertrand | A new characterization of threedimensional simple points[END_REF]. The algorithms used are detailed in Appendix A. Our topology preserving segmentation of the WM consists in performing a homotopic dilation of a seed set of voxels, called S, constrained to only add voxels from the set WM, knowing that S is topologically equivalent to a filled sphere. The result of this operation is denoted by SWM. For example, S could be made of single voxels chosen in the white matter, but we describe below a way to obtain a seed that is closer to the expected result, and thus leads to a more robust segmentation.

To obtain the seed S, we first compute a surface skeleton SK of WM, by dilating using Algorithm 2 as described in Appendix A. Then, we perform an homotopic erosion, constrained by SK, of a full cuboid that includes SK.

Finally, we perform an homotopic dilation of the same seed set S, constrained by the set SWM ∪ GM to only add GM and WM voxels, and we substract SWM from the result to obtain the corrected GM.

This method is performed on 3D sets, but for clarity we illustrate it on a 2D reduced example in Figure .2. Notice that small black components in

Partial volume relabelling and fractional content

The main contribution of the topology is the relabelling of missegmented GM voxels in hidden sulci as mixed GM/CSF. Once the topologically corrected WM, GM, CSF, WM/GM and GM/CSF segmentations are obtained, the portion of pure tissue, called here fractional content F , is computed for each mixed voxel by estimating the local contribution of each pure tissue.

We assume that each voxel contains at most two tissues and the new labelling corresponds only to the mixed voxels WM/GM and GM/CSF. For each mixed voxel, the fractional content F ranges between [0, 1] depending on the amount of pure tissue. Thus, for pure tissue voxels the fractional content F j are set to 1 for the class j and 0 otherwise. For mixed voxels (x ∈ W M/GM, GM/CSF ), the fractional content F j/k between both pure tissues j and k is computed using the intensity I(x) of the image and the robust local averages of the closest pure tissue voxels µ j (x) and µ k (x), such that:

F j/k (x) = U µ k (x) -I(x) µ k (x) -µ j (x) (1)
where U(•) is a limiter restricting the range of the fractional content to [0, 1].

Unlike [START_REF] Shattuck | Magnetic resonance image tissue classification using a partial volume model[END_REF] [START_REF] Coupe | An optimized blockwise nonlocal means denoising filter for 3-d magnetic resonance images[END_REF].

Pure tissue voxels are selected by eroding pure tissue segmentations using a 2mm radius, therefore reducing the influence of any mixed voxel. To illustrate the spatial differences in signal, GM intensity was measured over the population of 20 young adults scans, acquired as described in Subsection 3.3. Figure .6 shows the local average intensity of GM across the brain for an individual. In this example, precentral gyrus presented a higher average value than the temporal or occipital lobe. The same pattern appear in average in all the healthy individuals. The contrast between the tissues has been also measured using the Equation 2as

F Contrast = µ W M -µ GM µ GM -µ CSF (2)
where µ WM , µ GM and µ CSF are the regional averages of WM, GM and CSF respectively, which can be considered as a measure of the contrast between WM and GM normalized by the CSF. Figure .7 shows the regional differences for the population of 20 young controls.

INSERT FIGURE .7 HERE

Experiments

To evaluate our method, named hereafter as Topologically-corrected Partial Volume (TPV), we used different brain MR data sets including simulated and real images. The purpose was twofold, firstly to illustrate the effect of the topology correction in the estimation of fractional content for mixed voxels, and secondly to compare the obtained results with those publicly available in the area. After that, the method was integrated to our voxel-based cortical thickness estimation pipeline. Experiments demonstrated that the overall method showed a better estimate of thickness and a high reproducibility on real data.

Simulated MR data

A set of 15 simulated MR brain images was obtained from the BrainWeb who underwent 4 scans at baseline and 4 more scans during a subsequent session after a short delay (less than 90 days), were randomly selected from the Open Access Series of Imaging Studies (OASIS) database [START_REF] Marcus | Open access series of imaging studies (OASIS): Crosssectional MRI data in young, middle aged, nondemented, and demented older adults[END_REF], available at www.oasis-brains.org. For each session, an average motion-corrected image (co-registered average of all available data) was used for our tests. The scans were T1-weighted Magnetization Prepared RApid Gradient Echo (MP-RAGE) in sagittal orientation with isotropic 1mm 3 resolution (256 × 256 × 128 pixels). This data was used to assess the precision of the method when classifying partial volume voxels. We also tested the robustness when the method was integrated in our voxel-based cortical thickness estimation pipeline [START_REF] Acosta | Automated voxel-based 3D cortical thickness measurement in a combined lagrangian-eulerian PDE approach using partial volume maps[END_REF], particularly when the detection of deep sulci was improved.

Error and similarity measures

To quantitatively evaluate performance of the method over both simulated and real MR data sets and compare these results with other well-known results, we used two different metrics: the root mean square (RMS) error for comparison of PV classification maps, and the Jaccard similarity measure for comparison of the corresponding crisp tissue segmentations. The RMS error was used to quantify the differences between the fractional content calculated for each tissue and the corresponding values in the ground truth fuzzy membership images. As in [START_REF] Shattuck | Magnetic resonance image tissue classification using a partial volume model[END_REF], the RMS error between two images X and Y is calculated as

e RM S (X, Y ) = 1 |Ω| k∈Ω |y k -x k | 2
were Ω is the brain region, x k and y k are the image intensities at position k.

The Jaccard similarity metric, also known as the Tanimoto coefficient, measures the amount of overlap (agreement) between two images X and Y by taking the ratio between the size of their intersection and the size of their union:

J(X, Y ) = |X ∩ Y | |X ∪ Y |
This metric yields values between 0 and 1, where 0 means complete dissimilarity and 1 stands for identical images.

Results and discussion

BrainWeb

Performance of our TPV method was firstly assessed on the simulated brain images from BrainWeb. One example of the resulting PV maps for WM, GM and CSF, compared with the available ground truth, on the synthetic brain volume, 3% noise level and 20% bias field, is depicted in Figure .8.

Comparisons between our method and a classical MAP approach are shown in Figure .9 for the computed GMPVC fractional content map. It must be noted that compared to a classical MAP approach as in [START_REF] Shattuck | Magnetic resonance image tissue classification using a partial volume model[END_REF], the sulci were better delineated by introducing the topological constraints (Fig-

ure 9(g)). In this example, a deep sulci voxel with similar intensity to the average GM, will be classified as GM and not as a mixed GM/CSF voxel unless anatomical constraints are introduced. The mean RMS error of fractional content over the entire BrainWeb data set significantly decreased to 6.1% (p < 0.01) for the obtained GMPVC map, as compared with the results reported in [START_REF] Shattuck | Magnetic resonance image tissue classification using a partial volume model[END_REF] The variability between different regions in the brain may affect the performance of PV classifiers [START_REF] Chiverton | Adaptive partial volume classification of mri data[END_REF]. To illustrate this effect, we used the automated anatomical labeling (AAL) template [START_REF] Tzourio-Mazoyer | Automated anatomical labeling of activations in spm using a macroscopic anatomical parcellation of the mni mri single-subject brain[END_REF] to calculate the RMS error within each region as in [START_REF] Chiverton | Adaptive partial volume classification of mri data[END_REF]. Averaged results for different levels of noise are shown in Figure .10. As a low variability with respect to the bias field was observed, the depicted value corresponds to the average over all the bias field levels (0%, 20% and 40%). The smallest errors appeared in the amygdala (42xx), the insula (30xx), the supplementary motor area (24xx) and the olfactory (25xx); while lower agreement was found in the basal ganglia (70xx), the middle occipital (52xx) and the parietal superior (61xx).

INSERT FIGURE .10 HERE INSERT FIGURE .11 HERE

We also compared our TPV method with the results reported by Chiverton and Wells (2008) (GSR) and Shattuck et al. (2001) (SMAP). The results are depicted in Figure .11. Evidence suggests that the local average intensity strategy makes the classification more robust to bias field variations, and on average performs better than other methods for low levels of noise (1% to 7%) and bias field of 20%. We point out the fact that GSR does not explicitly take into account the bias field, hence its effect appears in the reported results.

Real MR Data

OASIS

The reproducibility was measured by applying the method to two of the MR scans from the same individual from the OASIS database. We compared the results with the MAP classifier as in [START_REF] Shattuck | Magnetic resonance image tissue classification using a partial volume model[END_REF]. Significant improvements in GM PV estimation were brought by the topology correction. The reproducibility error decreased by 8.8% in GM and 8.5% in WM (p < 0.001), measured as the RMS between the PV maps obtained on the rigidly registered baseline and repeat scans. Likewise, when comparing the crisp segmentations obtained by thresholding by 0.5 the baseline and repeat GM PV maps, the Jaccard similarity measure increased by 3.5% in GM. To compute crisp segmentations, each mixed voxel was assigned to the tissue class with the highest fractional content and the obtained segmentation were subsequently compared.

IBSR

Our method was also compared with both TMCD (trimmed minimum covariance determinant) [START_REF] Tohka | Fast and robust parameter estimation for statistical partial volume models in brain MRI[END_REF] and MMC (mixture model clus-tering) [START_REF] Noe | Partial volume segmentation of cerebral MRI scans with mixture model clustering[END_REF] on the IBSR data sets. Since the ground truth is available as manual segmentations performed by clinical experts, we compared the segmentations obtained from the crisped PV maps. Figure .12

shows an example of the ground truth provided by IBSR and a hard segmentation calculated after applying our method. were included for reference. Significant improvements in GM classification were demonstrated using the TPV, compared to a MAP classifier. The similarity measure (Jaccard) was improved by 8.7% in GM and 2.6% in WM (p < 0.001).

INSERT FIGURE .12 HERE

Poor similarity results were obtained in 5 cases, which exhibited strong shading artifacts that impeded a reliable GM and WM classification. Similar findings were presented in [START_REF] Noe | Partial volume segmentation of cerebral MRI scans with mixture model clustering[END_REF], who excluded them from the analysis. We also observed that the anisotropy in the images biased the computation of the local averages. Table .2 summarizes the mean (± standard deviation) of the Jaccard similarity values for each method, excluding the volumes with too severe intensity inhomogeneity. In average, our TPV method performed better for WM and GM compared to the others, excepting averaged GM segmentation against [START_REF] Noe | Partial volume segmentation of cerebral MRI scans with mixture model clustering[END_REF]. It must be noted that when the PV maps were used to generate the crisp segmentations, the mixed GM/CSF voxels in deep sulci with fractional content above 0.5 might be wrongly reclassified as GM. Under those conditions, the contribution of topology correction in the segmentation can not be fully and Thickness estimation with the proposed method (TPV) showed a higher reproducibility compared with the measure performed after partial volume classification using [START_REF] Shattuck | Magnetic resonance image tissue classification using a partial volume model[END_REF]. As can be seen in Figure .15, the differences in cortical thickness between scans were reduced after applying the TPV. The Pearson correlation coefficient was 0.915 in average and a paired ttest did not reveal any significant differences between the two measurements (p < 0.1). Also, the difference between scans was decreased by 13.7% in average, as shown in Table .3.

INSERT FIGURE .15 HERE INSERT TABLE .3 HERE

By using the proposed method, we found a mean (± std. dev.) cortical thickness over the whole brain of 2.08mm (± 0.11) for all the subjects, which is within the accepted range of cortical thickness for healthy young adults.

In previous studies, when the PV is not taken into account as in [START_REF] Yezzi | An Eulerian PDE approach for computing tissue thickness[END_REF], the computed mean thickness for the same population was 4.69mm (± 0.11). And when the PV classification method proposed by [START_REF] Shattuck | Magnetic resonance image tissue classification using a partial volume model[END_REF] is used, without any topology correction, the computed mean thickness was 3.06mm (± 0.25); using those same PV maps, but correcting the topology problems, decreases the mean thickness to 2.75mm (± 0.17). topology corrections and TPV, in average in this further step the number of voxels above 2.5mm has been reduced.

INSERT FIGURE .16 HERE

Conclusion

We (c) Differences between topology and TPV, in average the number of voxels above 2.5mm has been reduced consolidating the average thickness around 2.5 mm (typical value for young adults).

  presented a local adaptive Gradient-controlled spatial regularizer (GSR) using a Markov Random Field to model the class membership and a Markov chain Monte Carlo (MCMC) simulation to adapt the model to the observed data. The labelling error may remain high because the intensity inhomogeneities (not explicitly modelled)and the noise may lead to misdetection of mixed voxels mainly in tight sulci, representing a portion of GM/CSF/GM within the same voxel.

Figure 2

 2 Figure 2(b) can correspond to tunnels in the 3D image, thus simple connected component filtering would not give the correct region. Figures .3 and .4 show further examples in 3D.

  Finally, the computed averages are propagated back towards the location of the mixed voxels x, resulting in values of µ j (x) and µ k (x) that represent the average of the closest pure tissue voxels (Figure.5). The GM fractional content map is eventually defined as F GM/WM ∪ F GM ∪ F GM/CSF . Using a robust local mean overcomes issues related to intensity inhomogeneities and variations of pure tissue signal across the image, weighting accordingly the signal when computing the fractional content.INSERT FIGURE .5 HERE INSERT FIGURE .6 HERERegional differences in the cell structure and the distribution of different layers of the cortex result in variation of regional intensity differences for the same tissue across the brain. These differences produces local variation of contrast between the tissues that might be pronounced with ageing[START_REF] Salat | Age-associated alterations in cortical gray and white matter signal intensity and gray to white matter contrast[END_REF]. Whereas global homogeneity assumptions will bias the voxel fractional content estimation, a local computation of intensity averages for pure tissue yields a more accurate value, which accounts for the changes in cytoarchitecture visibles in MR. A local estimate allows also to overcome the issues of intensity inhomogeneities due to the artifacts during the acquisition.

  Figure 13(b) depicts the results of the comparison for the GM in the 20 normal subjects. As in (Chiverton and Wells, 2008), results of manual expert segmentation and pure tissue classsification presented by Ibrahim et al. (2006) (HMM, hidden Markov model)

Fig. . 16

 16 Fig. .16 depicts in histograms the impact of the topology correction and the accurate PV estimation on the cortical thickness calculation task. The higher thickness values produced after the first PV classification dissapeared when the topology of GM is corrected and the accurate PV value is computed

  have described a simple and fast technique to improve PV estimation of brain tissues from T1W MRI. It improves the detection of hidden mixed voxels in deep sulci by correcting for the topology errors in the segmentation and uses local averages to better estimate the fractional content. We show that fractional tissue content estimation can be improved for low levels of noise and regardless the intensity inhomogeneity, resulting in superior brain tissue segmentations. Topology correction improved the classification of mixed voxels in opposed banks of buried sulci by assuming GM as a continuous layer following the WM, with the topology of a filled sphere. Local modelling of tissue intensities helps to overcome the issues related with local intensity inhomogeneity and tissue MR properties across the image. Even with a preprocessing stage to correct the intensity inhomogeneities, pure cortical tissues show different intensity levels in the MRI. This suggests that the tissue properties are different depending on the region of the brain. Accuracy and precision were demonstrated and comparisons with other methods showed comparative performance with simulated and real MR data.We demonstrated the usefulness of the method to improve the accuracy of the cortical thickness estimation. By labelling mixed GM/CSF voxels in deep sulci and by recomputing a spatially compensated PV map, the measure of thickness in difficult regions is improved. Our method showed a high reproducibility on real data, with an extremely good agreement between the baseline and repeat scans. The computed values of thickness for young adults are similar to the ones reported previously in the literature. In the future, we plan to use our technique on clinical data to study cortical atrophy in Alzheimer's disease and other neurodegenerative diseases. We intend also to develop voxel-based techniques for inter-subject comparisons, a challenging issue given the large anatomical variability between patients.

  Figures and Legends586

Figure . 1 :Figure . 3 :

 13 Figure .1: Overall process for topology-corrected PV estimation in MR images

Figure . 5 :

 5 Figure .5: Schematic view of the local tissue averages for a given mixed voxel, where d1 and d2 relates to the closest voxels in the pure tissues.

Figure . 6 :

 6 Figure .6: Averaged intensity within the connected components of the pure GM, computed as the interquartile mean (IQM) within a 5mm radius sphere on an OASIS example data, normalized by the Maximum of intensity. The differences between the regions clearly appear. Thus, GM tissue intensity will be different between the regions and global homogeneity assumptions will slightly bias the computation of partial volume.

Figure . 7 :Figure . 8 :Figure . 9 :

 789 Figure .7: AAL template showing the regional differences incontrast between WM and GM over the surface, by calculating the ratio µWM-µGM µGM-µCSF . Darkest colours indicate bigger ratios, light colours indicate small values. Left: lateral and Right: medial views.

  Interslice distance is 3mm and the number of slices for each volume varies between 60 and 65. The data sets have various levels of artifacts,

	as low contrast and relatively large intensity gradients, that further affects
	performance of the algorithm. CMA also provides expert tissue labellings of
	each brain into WM, GM, and CSF, together with reference similarity values
	for some classification techniques.
	3.3. Cross sectional series of real MR scans
	20 young healthy subjects (12 female, 8 male; age between 19 -34 years),
	method on each combination of 1%, 3%, 5%, 7% or 9% noise levels together
	with 0%, 20% or 40% intensity nonuniformities. BrainWeb also provides the
	fuzzy tissue membership volumes, one for each tissue class, together with a
	discrete anatomical model of the simulated normal brain.
	3.2. Manually segmented real MR data
	20 normal MR brain data sets and their manual segmentations were ob-
	tained from the Internet Brain Segmentation Repository (IBSR), provided

Simulated Brain Database, maintained by the McConnell Brain Imaging Centre at the Montreal Neurological Institute

[START_REF] Cocosco | Brainweb: Online interface to a 3D MRI simulated brain database[END_REF]

and available at www.bic.mni.mcgill.ca/brainweb. Each simulation was a 1mm 3 isotropic T1-weighted MRI volume with dimensions 181 × 217 × 181, generated with varying noise level and intensity inhomogeneity. We tested our by the Center for Morphometric Analysis at Massachusetts General Hospital and available at www.cma.mgh.harvard.edu/ibsr. The data sets were acquired along the coronal axis with slice dimension of 256 × 256 and 1mm 2 resolution.

  . Overall, a good agreement was shown between the computed PV maps and the ground truth, available as fuzzy tissue membership volumes. RMS errors for different noise and intensity nonuniformity levels are shown in Table.1. As expected, the computed error was robust to the bias field, which additionally validates the local averaging approach rather than the global one.

	INSERT FIGURE .8 HERE
	INSERT FIGURE .9 HERE
	INSERT TABLE .1 HERE

accurately validated with this experiment. Nonetheless, we report these results for completeness.

INSERT TABLE .2 HERE INSERT FIGURE .13 HERE

Computational performance

On each image of the BrainWeb data set, after the initial MAP segmentation, the topology correction and PV fractional content estimation takes less than 10 minutes. For the OASIS data sets, the procedure takes about 9 minutes, while for the IBSR images the topology correction and PV fractional content estimation takes less than 4 minutes. Operations were encoded in a single-thread application and then executed in a standard Intel Core 2 Duo (3.00GHz, 2 GB RAM) machine running Linux.

Deep sulci cutting and cortical thickness estimation on real data

We integrated the proposed sulci detection and improved partial volume classification methods to our cortical thickness estimation pipeline [START_REF] Acosta | Automated voxel-based 3D cortical thickness measurement in a combined lagrangian-eulerian PDE approach using partial volume maps[END_REF], as depicted in Figure .14. Then, we computed the thickness, at two different acquisition times, for the same 20 young healthy subjects from the OASIS database [START_REF] Marcus | Open access series of imaging studies (OASIS): Crosssectional MRI data in young, middle aged, nondemented, and demented older adults[END_REF] Homotopic transformations are topology-preserving procedures that consist of sequentially deleting or adding simple points. This operation works only on binary images, such as the pure tissue segmentations, where each voxel is considered as a point. Informally, a simple point of an object X is a point that can be added or removed from X without changing the topological characteristics of X. It is possible to locally characterize simple points in 3D using two topological numbers T and T b [START_REF] Bertrand | A new characterization of threedimensional simple points[END_REF].

Thus, skipping some technical details, let A(x) be the set of points of X \ {x} lying in a neighborhood of x, and let Ab(x) be the set of points of the complement of X (background) lying in a neighborhood of x. Then, T (x) (resp. T b(x)) is the number of connected components of A(x) (resp. Ab(x)).

A point x is simple if and only if T (x) = T b(x) = 1. Topological numbers are useful for classifying points of an object X based on local topological characteristics: for example, a point x such that T b(x) > 1 characterizes a region of the object which separates (locally) its background into several parts.

Based on these notions, given an object X, a subset I of X and a priority function P , Algorithm 1 computes an homotopic erosion of X constrained by I, that is, an object that is topologically equivalent to X, that contains I and that has no simple point outside I. In this algorithm, the priority function P is usually chosen as the inverse of the distance to I, in order to select in the first place the points that are farthest to the set I. This choice will be assumed in the remaining operations.

Applying Algorithm 1 to the complementary sets of X and I, then in-