Damien Durville
email: damien.durville@ecp.fr

FINITE ELEMENT SIMULATION OF TEXTILE MATERIALS AT THE FIBER SCALE

A general approach to simulate the mechanical behavior of textile materials by taking into account all their constitutive elementary fibers and contacts between them is presented in this paper. A finite element code, based on an implicit solver, is develop to model samples of woven fabrics as assemblies of beams submitted to large displacements and developing contact-friction interactions between themselves. Special attention is paid to the detection and modeling of the numerous contacts occurring in this kind of structures. Robust models and efficient algorithms allow to consider samples made of few hundreds of fibers. The unknown initial configuration of the woven fabric is obtained by simulation according to the chosen weaving pattern, providing useful informations on fibers trajectories and tows cross-sections shapes. Various loading cases can then be applied to the textile sample, in presence of an elastic matrix, so as to identify its mechanical behavior under different solicitations. Results of simulations performed on coated fabric samples are given to illustrate the approach.

INTRODUCTION

Because the global behavior of textile composites depends on local phenomena occurring within the fabric at the scale of tows or fibers, a good understanding of mechanisms taking place at the scale of fibers is important. The nonlinear behavior of dry fabric or fabric-reinforced composites often originates in phenomena such as transverse compression of yarns or tows, sliding between tows or fibers, or locking between fibers. Predicting these phenomena is very hard without focusing the analysis at the scale of fibers.

A textile composite is made of components of different level. At a mesoscopic level, it can be seen as an assembly of yarns or tows, coupled with an elastic matrix. Modeling the structure at this mesoscopic level requires to be able identify the global behavior of tows or yarns. Because it involves nonlinear couplings between loads and strains in different directions due to interactions between fibers, the behavior of yarns is complex to identify. Going down to the scale of fibers allows to bypass the identification of intermediate models at upper levels and to base the global model only on mechanical properties of fibers and matrix and of contact-friction interactions between fibers.

Recent improvements in both computing capacity and modeling and algorithmic aspects, make the finite element method an appropriate tool to study textile composites. This paper presents a general approach, based on this method, to model the behavior of samples of dry fabric or textile composite materials, taking into account all elementary fibers together with the interactions developed between them and the coupling with an elastic matrix.

As an additional difficulty, when dealing with woven structures, the geometry of the initial configuration of the fabric (individual fiber trajectories) can not be known explicitly. Before starting with applying various loading cases on textile samples, the first task of the simulation is therefore to compute this initial geometry according with the weaving pattern.

The global problem is set as to find out an equilibrium solution for an assembly of fibers submitted to large displacements, developing contact-friction interactions, and bound to an elastic matrix. Fibers are represented by a 3D beam model, formulated in a large displacements and finite strains framework. The quasi-static incremental problem is solved using an implicit solver, which allows to consider large loading increments. The use of an implicit solver with many nonlinearities requires iterations until finding a solution satisfying the equilibrium and may lead to a high computation time. Because of the huge numbers of fibers and contacts, efficient and robust algorithms are needed to reduce this computation time. Some similar approaches can be found in the literature. At the scale of fibers, for the computation of the initial configuration for braided structures or 3D interlock woven fabric, digital elements have been used by Miao and al. [START_REF] Miao | Mechanics of textile composites: Micro-geometry[END_REF]. Since these digital elements have neither bending nor torsional stiffness, fibers must be tightened to find a solution. Finckh [START_REF] Finckh | Numerical simulation of mechanical properties of fabrics -weaving / Numerische Simulation der mechanischen Eigenschaften textiler Flächengebilde -Gewebeherstellung[END_REF] proposed to simulate the weaving process and to apply dynamic loading cases using an explicit solver. Other approaches tackle the problem at the scale of yarns, representing yarns by beams or 3D models and studying interactions between them [START_REF] Boisse | A mesoscopic approach for the simulation of woven fibre composite forming[END_REF][START_REF] Lomov | Meso-FE modelling of textile composites : Road-map, data flow and algorithms[END_REF][START_REF] Bridgens | Membrane material behaviour: concepts, practice and developments[END_REF]. The simulation approach developed in this paper has been previously applied to other kinds of entangled media [START_REF] Durville | Numerical simulation of entangled materials mechanical properties[END_REF][START_REF] Durville | Finite element simulation of the mechanical behaviour of textile composites at the mesoscopic scale of individual fibers[END_REF].

To introduce the global approach, a first section is dedicated to the beam model employed to represent fibers. The contact modeling, including geometrical and mechanical aspects, is then described. In a third section, the way the initial configuration is computed according to the weaving pattern is explained. After presenting how an automatically meshed elastic matrix can be added and coupled to the woven fabric, results of different loading tests on dry fabrics and composites, for both plain weave and twill weave are given.

3D BEAM MODEL

Introduction

Beam models are usually based on the assumption that their cross-sections remain rigid. The motions of these cross-section can therefore be described by two kinematical fields and six degrees of freedom : one translation of the center of the section, and one rotation. However, these usual models have two main disadvantages : the handling of large rotations requires complex formulations, and they can not account for deformation of cross-sections. To overcome these difficulties, we use a richer model described by nine degrees of freedom.

Enriched kinematical beam model

The enriched kinematical beam model is based on a Taylor expansion of the placement of particles of the beam with respect to the line of centroids. If we denote z (z 1 ,z 2 ,z 3 ,) a material particle defined in the reference configuration (see Fig. 1), the placement x (z) of this particle can be expressed as follows :

x (z 1 ,z 2 ,z 3 ) = x (0,0,z 3 ) + z 1 ∂ 1 x (0,0,z 3 ) + z 2 ∂ 2 x (0,0,z 3 ) + o(z 1 ,z 2 ).

[1]

Denoting respectively x 0 , g 1 and g 2 the three kinematical vectors introduced as first terms of the above expansion, we assume as kinematical model that the placement of any material particle of the beam is expressed as function of these vectors in the following way :

x (z 1 ,z 2 ,z 3 ) = x 0 (z 3 ) + z 1 g 1 (z 3 ) + z 2 g 2 (z 3 ). [START_REF] Finckh | Numerical simulation of mechanical properties of fabrics -weaving / Numerische Simulation der mechanischen Eigenschaften textiler Flächengebilde -Gewebeherstellung[END_REF] With this model, three vectors are used to describe any cross-section of the beam : one vector to give the position of its center, x 0 , and two other vectors, g 1 and g 2 , called section vectors, which determine the orientation of the section. Corresponding to this model for the placement, we assume that the displacement u (z) of any particle is also expressed by means of three vectors as follows :

u (z 1 ,z 2 ,z 3 ) = u 0 (z 3 ) + z 1 h 1 x (z 3 ) + z 2 h 2 (z 3 ). [3]
Figure 1. The three kinematical vectors used to define the placement in the kinematical beam model.

Since no conditions are set for the two section vectors, neither in norm, nor in angle, plane deformations of cross-section can be considered by the beam model, in addition to usual deformations (stretching, torsion, bending and shear).

Strains and constitutive law for the beam model

With the above kinematical beam model, because it is able to reproduce plane deformations of cross-section, no term of the Green-Lagrange strain tensor is constrained to be zero, and fully 3D effects, such as Poisson's effect (strains in transverse directions induced by stretching), can be caught. Because the strain tensor is fully three-dimensional, an usual 3D elastic constitutive law, taking into account strains in all directions, can be used to model the behavior of each beam.

CONTACT MODELING

Introduction

Contact modeling is crucial point since contact-friction interactions make the specificity of the behavior of textile structures. The main difficulties are related to the fact that contact areas can be numerous, and are constantly evolving, as contacts may appear or disappear anywhere and at any time.

Geometrical detection of contact

Since we are working in a large displacements framework, and because the rearrangement of fibers during the simulation, the process of contact detection needs to be regularly repeated during the calculation. For this reason, this process must be efficient and with low CPU cost.

Contact configurations between fibers can be very different depending on their relative positions. Very localized and reduced contact areas can be observed at crossings between fibers from two different tows, but there can also be long and continuous contact lines between two neighboring fibers of the same tow.

To face these different configuration, we suggest a general method to handle contact, based on the generation of contact elements, made of pairs of material particles attached to the surface of fibers, and based on geometrical criteria.

Whereas in usual methods, contact interactions are determined by fixing a first point on one of the structures and by searching a target point on the opposite geometry, using generally a normal direction to one surface, we suggest do define contact elements from an intermediate geometry.

This geometry has to be generated for each region where two parts of fibers are close enough and where contact is likely to produce. This way of searching contact from a third geometry offers a symmetrical treatment of contacting structures.

The first step to create contact elements is to determine in the global collection of beams, a set of proximity zones (see Fig. 2), defined as pairs of parts of lines holding beam meshes, which are close to each other. The determination of proximity zones is performed, for each couple of beams in the assembly, by distributing test points on one of the lines and calculating the nearest point on the opposite line.

Definition and generation of contact elements

We define a contact element E c (ij) (x G ) as a pair of material particles, z (i) and z (j) , belonging to two different beams (i) and (j), and predicted to enter into contact at a given location of the intermediate geometry :

E c (ij) (x G ) = (z (i) ,z (j) ) such that : x (z (i) ) = x (z (i) ) = x G . [4]
The equality between positions in the above equation must be understood in the sense of a prediction.

The discretization of the contact problem is performed by distributing on all intermediate geometries in the structure series of discrete locations where contact element are created. The fact that the construction of contact elements depends on the geometry of fibers, and therefore on the displacement solution of the problem, introduces in the formulation an additional nonlinearity, and requires the determination of contact elements to be repeated as the solution evolves.

Adding of a normal contact direction

In order to formulate linearized contact conditions, a normal direction, denoted N (E c (ij) ), is associated to each contact element E c (ij) . This direction sets the direction according to which distance between interacting particles will be measured.

Figure 6. Definition of a normal direction for a contact element.

This contact direction is generally taken as the one between the two centers of cross-sections involved in contact. Once again, such a definition introduces a nonlinearity in the global formulation since the normal contact direction depends on the solution itself.

Kinematical contact condition

In order to prevent interpenetration between beams, kinematical conditions are set for all contact elements. They aim at imposing that the distance between particles of a contact element, measured according to the normal contact direction, remains positive. This distance, defined as the gap, is expressed for each contact element as follows :

gap (E c (ij) ) = (x (z (i) ) -x (z (j) ),N (E c (ij) ) ≥ 0.

[5]

Constitutive laws for contact and friction

A constitutive law for normal contact is considered in order to link normal reactions to the gaps measured at contact elements. This constitutive law has the form of a penalty law, and is regularized by a quadratic part for very small penetrations in order to stabilize the contact algorithm (see Fig. 6). The penalty coefficient is adjusted for each proximity zone in order to limit the maximum penetration for each zone. As far as the friction is concerned, we use a Coulomb's law including a small reversible elastic displacement before sliding occurs.

Figure 7. Normal contact law and tangential friction law used for interactions.

Global algorithm to solve the nonlinear problem

To solve the problem for each loading increment, imbricated loops are introduced to iterate on the different nonlinearities. For each increment, with a first level loop, iterations are made on the process of contact determination to generate contact elements. Then, these contact elements being fixed, iterations are made at a second level on the normal contact directions. The third level loop is dedicated to iterations of a Newton-Raphson algorithm on all other nonlinearities (contact status, friction direction, finite strains). Because of these three embedded loops, the total number of iterations may be high if a good rate of convergence is not obtained for the Newton-Raphson algorithm. Models characteristics need to be carefully adjusted for this purpose.

COMPUTATION OF THE INITIAL CONFIGURATION

Way of determining the initial configuration

One important task of the simulation is to determine the unknown initial geometry of the woven structure. Instead of simulating step by step the actual weaving process, which would require a large amount of CPU time, the initial configuration is computed by starting from an arbitrary flat configuration where crossing tows are interpenetrating each other (Fig. 8), and moving gradually fibers until tows are moved apart from each other at crossings. For each crossing, the weaving pattern provides a superimposition order which indicates which tow should be above the other. The goal of the transient stage of computation of the initial configuration is to make this superimposition order fulfilled by fibers from different tows at crossings. To achieve this, while standard normal contact directions (see 3.4.1) are considered between fibers of the same tow, between fibers belonging to different tows, we take as normal contact direction a vertical direction oriented according to the superimposition order defined at the crossing. By this means fibers of crossing tows are moved step by step upon each other. Once tows do not interpenetrate any longer, standard normal contact directions are considered for all contact elements to find the equilibrium configuration.

Computation of the initial configuration for a plain weave and a twill weave sample

In the following, all applications are made for the same initial arrangement of 12 glass fiber tows (6 tows in both fill and warp directions), described on Table 1. Different tows are employed in the fill and warp direction so as to produce unbalanced fabrics. Two weaving patterns, corresponding respectively to a plain weave and a twill weave, are applied to this initial arrangement made of 408 fibers (Fig. 8). During the process of determination of the initial geometry, around 80 000 contact elements are generated.

Figure 8. Initial configuration of tows before weaving.

The view of slices at some steps during the process of determination of the initial configuration (Fig. 9) shows how fibers of different tows are gradually moved apart.

ADDITION OF AN ELASTIC MATRIX

In order to consider a composite sample, two layers of an elastic material are added on both sides of the fabric. The mesh of the layers is carried out automatically by the software. This mesh is relatively coarse compared to the size of the fibers, and overlaps slightly the external fibers of tows. Doing this way, the meshes of the fibers and of the matrix are non-conforming. To ensure the mechanical coupling between the fabric and the fibers, junction elements are generated in the overlapping region between the matrix and the fibers. These elements couple pairs of material particles, one belonging to a fiber and the other belonging to the matrix, by means of an elastic spring whose stiffness is calculated as function of the Young's modulus of the matrix.

Figure 12. View of the solid mesh generated for the matrix on both sides of the fabric.

SIMULATION OF VARIOUS LOADING CASES

Presentation

Various loading cases can be simulated by applying different loads (displacements or forces) on the borders of the reconstructed samples of dry fabric or textile composite. Rigid bodies are introduced in order to drive globally sets of fiber ends or tow ends. The nonlinear effect at the start of the loading, which is stronger for the twill weave sample, can be mainly explained by a change in the fabric geometry due to a deformation of tows and a possible rearrangement of fibers in tows. This can be observed on Fig. 14, where slices of the twill weave at the beginning and at the end of the test show the changes in tows cross-section shapes and tows trajectories.

CONCLUSIONS

A general approach for the simulation by the finite element method of the mechanical behavior of fabrics and fabric-reinforced composites has been presented. This approach, formulated in a large displacements and finite strains framework, is based on the representation of all fibers constituting woven structures by means of 3D beam models, and on the taking into account of contact-friction interactions between fibers.

The detection of the numerous contacts occurring in general assemblies of fibers is one of the key points of the model. Based on the determination at a first level of proximity zones between fibers, and on the construction of intermediate geometries to approximate the actual contact zones, the method generates automatically contact elements made of pairs of material particles. This method is general enough to apply to many contact configurations as encountered between fibers in woven materials. Since the process of determination of contact elements depends on the relative positions of fibers, it has to be repeated during the solving of the problem for each loading step.

Robust models and efficient algorithms are required to get a good convergence rate for the solving of the nonlinear problem. Thanks to these optimized models and algorithms, samples of woven fabrics involving few hundreds of fibers and about 100 000 contact elements can be studied with a reasonable CPU time.

The model is first applied to the determination of the unknown initial geometry, by making the superimposition order at crossings defined by the weaving pattern be gradually fulfilled by fibers. Meaningful informations related to the yarns trajectories and shapes, and to the trajectories and curvatures of fibers, are obtained as results of this initial simulation.

In order to create a sample of textile composite, an elastic matrix is added to the numerically manufactured fabric, and is automatically meshed. Various loading cases can be than applied to the textile composite sample in order to identify its behavior under different kind of solicitations. Thanks to the nonlinearities considered in the model, and its ability to handle large loading increments, typical nonlinear effects of the complex behavior of textile materials can be rendered for a wide range of loadings.

The proposed simulation, based on the identification of very few parameters (mechanical properties of fibers, definition of tows arrangement and of the weaving pattern), appears as a suited tool to describe and to understand complex phenomena occurring at the scale of fibers, and to predict the complex global mechanical behavior of textile composite materials at a macroscopic scale, as well as very localized phenomena such as breakings at the microscopic scale of fibers.

Figure 2 .

 2 Figure 2. Definition of intermediate geometries based on proximity zones.The intermediate geometry is then defined, for each proximity zone, as the average between the two close parts of lines. Contact elements are generated from this intermediate geometry, using planes orthogonal to this geometry to determine couples of beam cross-sections candidates to contact. Material particles are then accurately located on the border of these cross-sections.

Figure 3 .

 3 Figure 3. Determination of couples of cross-sections candidates to contact.

Figure 4 .

 4 Figure 4. Selection of material particles on the borders of cross-sections to create a contact element.

Figure 5 .

 5 Figure 5. Contact elements (red spots) generated in a woven structure 3.4 Formulation of contact conditions

Figure 9 .

 9 Figure 9. Slices at different steps during the process of determination of the initial configuration for the plain weave and the twill weave. Many informations about the fabrics geometry are obtained as results of this initial simulation. Trajectories of both fibers and tows, varying shapes of tows cross-sections and local curvatures along fibers are very useful data which can be recovered from this first computation.

Figure 10 .

 10 Figure 10. Computed initial configurations computed for plain and twill weaves.

Figure 11 .

 11 Figure 11. Details of the computed initial configurations for plain and twill weaves.

6. 2

 2 Biaxial tensile tests on dry fabric samplesBiaxial tensile tests are applied on the plain weave and on the twill weave samples. For these tests, two different ratios have been considered for elongations in warp and fill directions : a ratio of 1, with equal elongations in both directions, and a ratio of 0, with no elongation in warp direction. The corresponding loading curves are plotted on Fig.13. They show the typical nonlinear behavior characterizing this kind of material.

Figure 13 .

 13 Figure 13. Loading curves in fill and warp directions for biaxial tensile tests on samples of unbalanced dry fabric for plain weave and twill weave.

Figure 14 .

 14 Figure 14. Slices of the twill weave fabric before (above) and after (below) a 2 % elongation test,showing deformation and rearrangement within tows cross-sections.6.3 Shear test on the twill fabric compositeA shear test is performed on the composite textile sample for twill weave. The view of the deformed mesh is presented on Fig.15. The shear deformation is possible until reaching a kind of locking when neighboring parallel tows come to contact. On a zoom view on Fig.16, the deformation of tows can be observed.

Figure 15 .

 15 Figure 15. View of the textile composite sample submitted to a shear test.

Figure 16 .

 16 Figure 16. Slices at the beginning and at the end of the shear test showing changes is the shape of tows cross-sections. 6.4 Bending test on the twill fabric composite The last test conducted on the composite sample is a bending test, for which opposite rotations are imposed to two borders of the sample. A global view of the sample is presented on Fig. 17, and a corresponding slice on Fig. 18. This test shows the ability of the model to handle large displacements.

Figure 17 .

 17 Figure 17. Global view of the composite sample after bending.

Figure 17 .

 17 Figure 17. Slice of the composite sample submitted to bending.

Table 1 .

 1 Characteristics of the studied samples.

	Characteristics of the studied samples
	Number of fibers in fill tows	44
	Number of fibers in warp tows	22
	Total number of fibers	408
	Number of nodes	35 112
	Number of dofs	316 008
	Number of contact elements	around 80 000


S Conforti

Di Lorenzo

J E Campagne

F Dulucq

C De

La Taille

G Martin-Chassard

M El

PARISROC is a complete read out chip, in AMS SiGe 0.35 m technology, for photomultipliers array. It allows triggerless acquisition for next generation neutrino experiments and it belongs to an R&D program funded by the French national agency for research (ANR) called PMm 2 : "Innovative electronics for photodetectors array used in High Energy Physics and Astroparticles" (ref.ANR-06-BLAN-0186). The ASIC (Application Specific Integrated Circuit) integrates 16 independent and auto triggered channels with variable gain and provides charge and time measurement by a Wilkinson ADC (Analog to Digital Converter) and a 24-bit Counter. The charge measurement should be performed from 1 up to 300 photoelectrons (p.e.) with a good linearity. The time measurement allowed to a coarse time with a 24-bit counter at 10 MHz and a fine time on a 100ns ramp to achieve a resolution of 1 ns. The ASIC sends out only the relevant data through network cables to the central data storage. This paper describes the front-end electronics ASIC called PARISROC.

I. INTRODUCTION

T he next generation of proton decay and neutrino experiments, the post-SuperKamiokande detectors as those that will take place in megaton size water tanks, will require very large surfaces of photodetection, even with large hemispherical photomultiplier tubes (PMT), and a large volume of data: the expected number of channels should reach hundreds of thousands. The PMm 2 project [1] proposes to segment the large surface of photodetection in macro pixel consisting of an array (2 m × 2m) of 16 photomultipliers (16 PMTs as SuperKamiokande module structure) connected to an autonomous front-end electronics (Fig. 1) and powered by a common High Voltage (HV). This R&D involves three French laboratories (LAL Orsay, LAPP Annecy, IPN Orsay) and ULB Bruxells for the DAQ. It is funded by the French National Agency for Research (ANR) under the reference ANR-06-BLAN-0186.

The micro-electronics group's (OMEGA from the LAL at Orsay) is in charge of the design and tests of the readout chip *Omega (Orsay Microelectronics Group Associated)/LAL/IN2P3 (http://omega.in2p3.fr/). Manuscript received November 13, 2009.

named PARISROC which stands for Photomultiplier ARrray Integrated in Si-Ge Read Out Chip.

The detectors are large tanks covered by a significant number of large photomultipliers (20-inch); the next generation neutrino experiments will require a bigger surface of photo detection and thus more photomultipliers. As a consequence the total cost has an important relief [START_REF] Genolini | PMm 2 : large photomultipliers and innovative electronics for next generation neutrino experiments[END_REF]. The project proposes to reduce the costs using:

1. 12-inch PMts, instead of 20-inch, with an improved cost per unit of surface area and detected p.e (due to a different industrial fabrication of the PMTs); 2. A smaller number of electronics, thanks to the 16 PMTs macropixel with a common electronics; 3. The designing of the ASIC as a System on a Chip (SoC) that processing the analog signal up to digitization (cost improvement only in a large scale of production); 4. A common High Voltage for the 16 PMTs so a reduced number of underwater cables, cables that are also used to brought the DATA to the surface; 5. The front-end closed to the PMTs that allow a suppression of underwater connectors.

PARISROC, a Photomultiplier Array Integrated

Read Out Chip.

The general principle of PMm 2 project is that the ASIC and a FPGA (Field-Programmable gate array) manage the dialog between the PMTs and the surface controller (Fig. 2).

Through the 100 m of data cable (used also to transmit the slow control, the power supply, the 10 MHz reference clock) the digitized data are transmitted with a dedicated protocol at a rate of 5 Mbit/s. The data are processed at the surface level on a network: the interface with the front end is carried out by a surface controller [START_REF] Genolini | PMm 2 : large photomultipliers and innovative electronics for next generation neutrino experiments[END_REF].

Alternative options may be chosen considering an analysis of the risks of this full underwater strategy. One option is that the Front-End electronics can be used in a traditional schema with the electronic "in surface". PARISROC can be integrated perfectly in a surface scheme.

II. REQUIREMENTS

The physics events researched in the detectors, as SuperKamiokande one, produce Cerenkov light that is spread over the PMTs. This kind of experiences we are looking for rare events and the number of pe per deposit MeV in the water is of the order of 10pe/MeV and may be viewed by O(10000) PMTs. So, for few MeV events, the small number of p.e is spread over a large number of PMTs and as consequence it is necessary being fully efficient to detect a single photo-electron (p.e). For large energy events (atmospheric or beam neutrinos) it is shown by simulation that the dynamic range for a single PMT should cover up to few hundred p.e (~ 300pe). This induces an electronic that covers a dynamic range from 1 to 300 pe and a trigger at 1/3 p.e.

All the type of events considered must be registered without any direct external trigger, this is called 'triggerless mode' and requires an electronic with autotrigger channels.

A precise time stamp of each event is required to reconstruct the topology of the events and so to synchronize the events among PMTs in each array and among the different arrays. This aspect brought to a requirement: an electronic with full independent channels. The most demanding in term of timing is the vertex reconstruction that needs typically 1 ns resolution considering the single electron jitter around 3 ns (for 12-inch PMT).

The mean PMT gain is about 3×10 6 but a spread is foreseen among the PMTs of the same array as a single HV is used. It was estimate from the production of the Auger experiment PMTs that the gain dispersion, at a given voltage is such that, the ratio between the highest and the lowest gain, is not more than 12 [START_REF] Genolini | PMm 2 : large photomultipliers and innovative electronics for next generation neutrino experiments[END_REF]. It is possible for the manufacturer to sort the PMTs, at a reasonable cost, when they are produced at a very large scale, with a gain ratio reduced to 6 in a batch of 16 PMTs. To compensate this not homogeneity a preamplifier with a variable and adjustable gain is required.

To summarize, the electronic requirements are: 1pe full efficiency, triggerless acquisition, 1ns of time resolution, 300 pe of dynamic range, high granularity, scalability, low cost, independent channels, water-tight, common high voltage, only digitized data to be send to the surface (DATA + VCC).

III. PARISROC ARCHITECTURE

A. Global architecture. The ASIC PARISROC (Fig. 3) is composed of 16 analog channels managed by a common digital part. Each analog channel (Fig. 4) is made of a voltage preamplifier with variable and adjustable gain. The variable gain is common for all channels and it can change, thanks to the input variable capacitance, on 3 bits. The gain is also tuneable channel by channel to adjust the input PMTs gain non-homogeneity, thanks to the switched feedback capacitance, on 8 bits. The preamplifier is followed by a slow channel, for the charge measurement, in parallel with a fast channel, for the trigger output. The slow channel is made by a slow shaper followed by an analogue memory, with a depth of 2, to provide a linear charge measurement up to 50 pC; this charge is converted by a Wilkinson ADC (8, 10, 12-bit). One follower OTA is added to deliver an analogue multiplexed charge measurement. The fast channel consists in a fast shaper (15 ns) followed by 2 low offset discriminators, to auto-trig down to 50 fC. The thresholds are loaded by 2 internal 10-bit DACs (Digital to Analog Converter) common for the 16 channels and an individual 4bit DAC for one discriminator. The 2 discriminator outputs are multiplexed to provide only 16 trigger outputs. Each output trigger is latched to hold the state of the response until the end of the clock cycle. It is also delayed to open the hold switch at the maximum of the slow shaper. An "OR" of the 16 trigger gives a 17th output. For each channel, a fine time measurement is made by an analogue memory, with depth of 2, which samples a 12-bit ramp, common for all channels, at the same time of the charge. This time is then converted by a Wilkinson ADC. The two ADC discriminators have a common ramp, of 8, 10, 12 bits, as threshold to convert the charge and the fine time. In addition a bandgap bloc provides all voltage references.

B. Digital part.

The digital part of PARISROC is built around 4 modules which are acquisition, conversion, readout and top manager [START_REF] Dulucq | Digital part of PARISROC: a photomultiplier array readout chip[END_REF]. Actually, PARISROC is based on 2 memories. During acquisition, discriminated analog signals are stored into an analog memory (the SCA: switched capacitor array). The analog to digital conversion module converts analog charges and times from SCA into 12 bits digital values. These digital values are saved into registers (RAM). At the end of the cycle, the RAM is readout by an external system. The block diagram is given below on Fig. 5. The first module in the sequence is the acquisition which is dedicated to charge and fine time measurements. It manages the SCA where charge and fine time are stored as a voltage like. It also integrates the coarse time measurement thanks to a 24-bit gray counter with a resolution of 100 ns. Each channel has a depth of 2 for the SCA and they are managed individually. Besides, SCA is treated like a FIFO memory: analog voltage can be written, read and erased from this memory.

Then, the conversion module converts analog values stored in the SCA (charge and fine time) in digital ones thanks to a 12-bit Wilkinson ADC. The counter clock frequency is 40 MHz, it implies a maximum ADC conversion time of 103 µs when it overflows. This module makes 32 conversions in 1 run (16 charges and 16 fine times).

Finally, the readout module permits to empty all the registers to an external system. As it will only transfer hit channels, this module tag each frame with its channel number: it works as a selective readout. The pattern used is composed of 4 data: 4-bit channel number, 24-bit coarse time, 12-bit charge and 12-bit fine time. The total length of one frame is 52 bits. The maximum readout time appears when all channels are hit. About 832 bits of data are transferred to the concentrator with a 10 MHz clock: the readout takes about 100 µs with 1 µs between 2 frames.

C. PARISROC layout

Figure 6 presents the PARISROC layout [START_REF] Martin-Chassard | PMm 2 ASIC: PARISROC, TWEPP 08 conference[END_REF]. The circuit has been designed in AMS SiGe 0.35 m technology [5]. The die has a surface of 17 mm 2 (5 mm X 3.4 mm) and will be package in CQFP160 case.

IV. MEASUREMENTS AND SIMULATIONS.

A. General tests.

A dedicated test board has been designed and realized for testing the ASIC (Fig. 7). Its aim is to allow the characterization of the chip and the communication between photomultipliers and ASIC. This is possible thanks to a dedicated Labview program that allows to send the ASIC configuration (slow control parameters, ASIC parameters, etc) and to receive the output bits via an USB cable connected to the test board. The Labview is developed by the LAL "Tests group". A signal generator is used to create the input charge injected in the ASIC. The signal injected is similar, as possible, to the PMT signal. In Fig. 8 is represented the input signal and its characteristics. The input signal, used in measurements and simulation, is a triangle signal with 5 ns rise and fall time and 5 ns of duration. This current signal is sent to an external resistor (50 Ohms) and varies from 0 to 5 mA in order to simulate a PMT charge from 0 to 50 pC which represents 0 to 300 p.e. when the PM gain is 10 6 . There is a good agreement between measurements and simulations in analog part results except for the noise values. To characterize the noise, the Signal to Noise Ratio (SNR) is calculated with reference to the MIP (1 p.e.). The noise differences are immediately evident: an additional low frequency noise is present in measurement (is now under investigation even if it is supposed to be tied to the power supply noise). A small difference has been noticed in measurement without the USB cable that allowed the communication between the test board and the Labwiev program: an rms noise value of 660 V (0.132 p.e.) for preamplifier and so a SNR value of 8 (slightly nearer to the simulation value of 12).

Another important characteristic is the linearity. The preamplifier linearity in function of variable feedback capacitor value with an input charge of 10 p.e. and with residuals from -2.5 to 1.35 % is represented on Fig. 9. The gain adjustment linearity is good at 2% on 8 bits. In order to investigate the homogeneity among the whole chip, essential for a multichannel ASIC, for different preamplifier gains are plotted the maximum voltage values for all channels. On Fig. 11 is given the gain uniformity. A good dispersion of 0.5%, 1.4% and 1.2% have respectively been obtained for gain 8, 4 and 2. This represents a goal for the ASIC.

B. DAC Linearity

The DAC linearity (Fig. 12) has been measured and it consists in measuring the voltage DAC (Vdac) amplitude obtained for different DAC register values. Figure 12 gives the evolution of Vdac as a function of the register for the DAC and the residuals with values from -0.1% to 0.1%.

C. Trigger Output.

The trigger output behaviour was studied scanning the threshold for different injected charges. At first no charge was injected which corresponds to measure the fast shaper pedestal. The result is represented on Fig. 13 for each channel. The 16 curves (called s-curves because of their shape) are superimposed that meaning good homogeneity. The spread is of one DAC count (LSB DAC=1.78 mV) equivalent to 0.06 p.e. The trigger output is studied also by scanning the threshold for a fixed channel and changing the injected charge. Figure 15 shows the trigger efficiency versus the DAC unit with an injected charge from 0 to 300 p.e. and on Fig. 16 is plotted the threshold versus the injected charge but only until 3 pe. In Fig. 17 are plotted the 50% trigger efficiency values, extracted from the plot in Fig. 16, converted in mV versus the injected charges. A noise of 10 fC has been extrapolated. Therefore the threshold is limited to 10  noise due to the discriminator coupling.

D. ADC.

The ADC performance has been studied alone and with the whole chain. Injecting to the ADC input directly DC voltages by the internal DAC (in order to have a voltage level as stable as possible) the ADC values for all channels have been measured. The measurement is repeated 10000 times for each channel and in the first panel of the Labview window (Fig. 18) the minimal, maximal and mean values, over all acquisitions, for each channel are plotted. In the second panel there is the rms charge value versus channel number with a value in the range 0.5 to 1 ADC unit. Finally the third panel shows an example of charge amplitude distribution for a single channel; a spread of 5 ADC counts is obtained. The ADC is suited to a multichannel conversion so the uniformity and linearity are studied in order to characterize the ADC behaviour. On Fig. 19 is represented the ADC transfer function for the 10-bit ADC versus the input voltage level. All channels are represented and have plots superimposed. This plot shows the good ADC uniformity among the 16 channels. In Fig. 20 is shown the 12-bit ADC linearity plots with the 25 measurements made at each input voltage level. The average ADC count value is plotted versus the input signal. The residuals from -1.5 to 0.9 ADC units for the 12-bit ADC; from -0.5 to 0.4 for the 10-bit ADC and from -0.5 to 0.5 for the 8-bit ADC prove the good ADC behaviour in terms of Integral non linearity. Once the ADC performances have been tested separately, the measurements are performed on the complete chain. The results of the input signal auto triggered, held in the T&H and converted in the ADC are illustrated in Fig. 21 where are plotted the 10-bit ADC counts in function of the variable input charge (up to 50 p.e). A nice linearity of 1.4% and a noise of 6 ADC units are obtained.

Fig . 1 .

 1 Fig .1. Principal of PMm 2 proposal for megaton scale Cerenkov water tank with 16 PMTs array shown in the bottom side of the picture and are emphasized the surface controller and the 100m of cable.

Fig. 2 .

 2 Fig. 2. Principle of the PMm 2 project. The main blocks are shown: the integration board that receive the PMT signal and transmit by 100 m of cable the 52 bits of DATA (for each channel) to the surface controller.

Fig. 3 .

 3 Fig. 3. PARISROC global schematic.

Fig. 4 .

 4 Fig. 4. PARISROC one channel analog part schematic. In the picture are emphasized, on the left the preamplifier with variable input and feedback capacitances; on the top, the Track &Hold block with the two capacitances to memorize charge and time measurements and the two ADC, with a common ADC ramp, to convert the analog values.

Fig. 5 .

 5 Fig. 5. Block diagram of the digital part.The sequence is made thanks to the top manager module which controls the 3 other ones. When 1 or more channels are hit, it starts ADC conversion and then the readout of digitized data. The maximum cycle length is about 200 µs. During conversion and readout, acquisition is never stopped. It means that discriminated analog signals can be stored in the SCA at any time.The first module in the sequence is the acquisition which is dedicated to charge and fine time measurements. It manages the SCA where charge and fine time are stored as a voltage like. It also integrates the coarse time measurement thanks to a 24-bit gray counter with a resolution of 100 ns. Each channel has a depth of 2 for the SCA and they are managed individually. Besides, SCA is treated like a FIFO memory: analog voltage can be written, read and erased from this memory.Then, the conversion module converts analog values stored in the SCA (charge and fine time) in digital ones thanks to a 12-bit Wilkinson ADC. The counter clock frequency is 40 MHz, it implies a maximum ADC conversion time of 103 µs when it overflows. This module makes 32 conversions in 1 run (16 charges and 16 fine times).Finally, the readout module permits to empty all the registers to an external system. As it will only transfer hit channels, this module tag each frame with its channel number: it works as a selective readout. The pattern used is composed of 4 data: 4-bit channel number, 24-bit coarse time, 12-bit charge and 12-bit fine time. The total length of one frame is 52 bits. The maximum readout time appears when all channels are hit. About 832 bits of data are transferred to the concentrator with a 10 MHz clock: the readout takes about 100 µs with 1 µs between 2 frames.

Fig. 6 .

 6 Fig. 6. PARISROC layout. The PARISROC has been submitted in June 2008; a first batch of 6 ASICs has been produced and received in December 2008 (a second batch of 14 ASICs in May 2009).

Fig. 7 .

 7 Fig. 7. Test Board.

Fig. 8 .

 8 Fig. 8. Input signal used for measurements and simulations.

Fig. 9 .

 9 Fig. 9. Preamplifier linearity vs feedback capacitor value.

Figure 10

 10 Figure 10 represents the slow shaper linearity for a time constant of 50 ns and a preamplifier gain of 8. The slow shaper output voltage in function of the input injected charge is plotted. Good linearity performances are obtained with residuals better than ± 1%.

Fig. 10 .

 10 Fig. 10. Slow shaper linearity; =50 ns and Gpa=8.

Fig. 11 .

 11 Fig. 11. Gain uniformity for Gpa = 8, 4, 2.



  Fig. 12. DAC Linearity.

Fig. 13 .

 13 Fig. 13. Pedestal S-Curves for channel 1 to 16.

Fig. 14

 14 Fig. 14. S-Curve for input of 10 p.e. for channel 1 to 16.

Fig. 15 .

 15 Fig. 15. Trigger efficiency vs DAC count up to 300 p.e.

Fig. 16 .

 16 Fig. 16. Trigger efficiency vs DAC count until 3 p.e.

Fig. 17 .

 17 Fig. 17. Threshold vs injected charge until 500 fC.

Fig. 18 .

 18 Fig. 18. ADC measurements with DC input 1.45V.

Fig. 19 .

 19 Fig. 19. 10-bit ADC transfer function vs input charge.

Fig. 20

 20 Fig. 20. 12-bit ADC linearity.

Fig. 21 .

 21 Fig. 21. 10-bit ADC linearity.

Table I

 I lists the simulation and measurements results for the three main blocks of the analog part: Preamplifier, Slow shaper and Fast shaper.

TABLE I .

 I ANALOG PART MEASUREMENTS

		Preamplifier Slow Shaper	Fast Shaper
		Gain P.A.=8 RC=50ns
		Meas./Sim. Meas./Sim.	Meas./Sim.
	Voltage	5mV/5.43mV 12mV/19mV	30mV/39mV
	( 1 p.e.)
	Rms noise	1mV/468V	4mV/2.3mV	2.5mV/2.4mV
	Noise in p.e. 0.2/0.086	0.3/0.125	0.08/0.06
	SNR	5/12	3/8	12/16


Maxim Fradkin

Cybèle Ciofolo

Benoît Mory

Gilion Hautvast

Marcel Breeuwer

Fully automatic segmentation of short and long axis cine cardiac MR

Introduction

Quantitative analysis of cardiac function requires delineation of the left ventricle (LV) in cine cardiac MR (CMR). Typically, this is done using short-axis (SA) images, however, acquisition of several long-axis (LA) views has become quite common. The latter can be used for the accurate and reproducible determination of the basal SA slice, known as one of the major inter-observer variability factors in SA LV measurements [1,2]. Since manual LV delineation is very tedious and time-consuming, automatic segmentation methods, enabling to obtain reproducible LV measurements, are highly desirable. Though many publications proposed such methods, mainly for SA, only very few of them reported acceptable accuracy.

Purpose

We propose a fully automatic method for delineation of the endo-and epicardial contours in SA and LA cine CMR images in order to provide automatic, accurate quantitative leftventricular functional assessment.

Methods

Our segmentation workflow consists of two steps: (1) automatic detection and delineation of the LV myocardium and valve plane in LA image(s); (2) segmentation of SA myocardial contours, using the results of step 1. For modelling the LV, we use ribbon-like template, relying on interpolating splines and defined by a very few nodes; its shape is controlled by nodes' position and the template width at each node (Fig. 1). The template is deformed using greedy optimization framework [3]. The optimization criterion comprises three terms: shape, contrast and region. The shape term controls the smoothness of the template as well as its resemblance to a pre-defined shape. The contrast term attracts the template boundaries to the myocardial borders in the image. Finally, the region term imposes homogeneity of the myocardium and the blood-pool regions in the underlying image. The optimization is done in coarse-to-fine manner, by relaxing the deformation type (from rigid to local) and adjusting the influence (weights) of different terms. The deformation framework is exactly the same for both LA and SA, the only difference is in the initial template positioning and deformation parameters. The initial position of the LA template is derived from the image acquisition geometry, while that for the SA templates is obtained by projecting the LA segmentation contours on the corresponding SA slices. We quantitatively assessed the performance of the method on a database of 35 cine CMR acquisitions (Philips Gyroscan NT Intera 1.5T, SSFP protocol, TE=1.6ms, TR=3.3ms, flip angle 50 0 for SA and TE=1.5-1.7ms, TR=3.0-3.4ms, flip angle 55 0 for LA). The SA ones included 9 to 14 slices (325 slices in total), and the LA ones 1 to 2 views (40 views in total). The assessment was made by comparing the automatically segmented contours with Gold Standard manual delineations provided by an expert.

Results

For SA images, the mean positioning error (MPS) was 1.3 ± 0.5mm for the endocardium and 1.5 ± 0.8mm for the epicardium, for a pixel size of 1.4-1.8mm. For LA images, the MPS was 1.3 ± 0.4mm (endocardium) and 1.1 ± 0.4mm (epicardium), for a pixel size of 1.3-1.8mm. Both results are under one pixel accuracy and comparable with inter-and intra-observer variability. An excellent correlation was found between "manual" and "automatic" LV volumes (r=0.98, P < 0.001), with mean difference -2.5 ± 7.6 ml (5.5 ± 3.9%). Note that all results were obtained without any user interaction. An example of segmentation results for one patient (8 SA slices and 2 LA views) is shown in Fig. 2. We validated the method for the end diastolic (ED) phase only, since the segmentation of the whole cardiac cycle can be then obtained using automatic propagation of ED contours [4]. In earlier investigations, we already showed that propagation preserves the segmentation accuracy within acceptable ranges while being much faster than phase-by-phase segmentation.

Conclusion

We presented a robust, accurate and efficient method for the fully automatic delineation of the myocardium contours in LA + SA cine CMR images, which can be used for accurate LV functional assessment.

[1] Weaver, A., et al.: Journal of Cardiovascular Magnetic Resonance 1, pp. 82-83 (2006) [2] Anderson, J., et al.:., The Open General and Internal Medicine Journal 1, 6-12 (2007) [3] Williams, D.J., Shah, M., CVGIP: Image Understanding 55(1), 14-26 (1992) [4] Hautvast et. al., IEEE. Trans. Med. Imag. vol. 25(11), pp. 1472-1482, 2006 Figure 1: Deformable templates modeling LA (left) and SA (right): centerline (dashed) is defined by a few nodes, the template width (grey ribbon) may vary from node to node.

Figure 2 :

 2 Figure 2: LV segmentation example for one patient: (from left to right, top to bottom): consecutive slices of SA volume from valve to apex; LA 4 chambers and 2 chambers.


A variational model for image texture identification R. Echegut and L. Piffet Université d'Orléans, Département de Mathématiques (MAPMO), UFR Sciences, Bâtiment de mathématiques, Route de Chartres B.P. 6759, 45067 Orléans cedex 2, FRANCE echegut.romain@wanadoo.fr, Loic.Piffet@univ-orleans.fr Variational models in image processings have been extensively studied during the past decade. There are used for segmentation processes (geodesic or geometric contours) and restoration and textures extraction purpose as well. Roughly speaking image restoration problems are severely ill posed and a Tychonov-like regularization process is needed. The general form of such models consists in the mimization of a "energy" functional :

F(u) = u -u d X + R(u) , u ∈ Y ⊂ X ,
where X, Y are (real) Banach spaces, R is a regularization operator, u d is the observed (or measured) image and u is the image to recover. The first term is the fitting data term and the second one permits to get a problem which is no longer ill posed. The most famous model is the Rudin-Osher-Fatemi denoising model (see [START_REF] Osher | Nonlinear total variation based noise removal algorithms[END_REF]) which leads to the following problem:

min u∈BV (Ω) λ|u| BV (Ω) + 1 2 u d -u 2 L 2 (Ω) . (P ROF )
Here, we are interested in textures extraction and/or image denoising. Recently people considered that an image can be decomposed in many components, each component describing a particular property . of the image (see [START_REF] Aubert | Image decomposition into a bounded variation component and an oscillating component[END_REF][START_REF] Osher | Image denoising and decomposition with total variation minimization and oscillatory functions. Special issue on mathematics and image analysis[END_REF][START_REF] Yin | A comparison of three total variation based texture extraction models[END_REF] for example). We agree with this idea and assume that the image f we want to recover from the data u d can be decomposed as f = u + v where u and v are functions that belong to different functional spaces: u is the "texture" part which involves (periodic or not) details (and noise as well) while v is a more regular part (usually called the "cartoon" component).

In a first section, we present the functional framework, introducing the BV 2 space, and the general variational model we consider. In section 2, we focus on numerical implementation and present the discretization process. Numerical tests are reported in the last section.

1 Functional framework and model

The BV 2 (Ω) space

Let Ω be an open bounded subset of R n , n ≥ 2 (practically n = 2). Following Demengel [START_REF] Demengel | Fonctions à hessien borné[END_REF], we define the space of bounded hessian functions that we call BV 2 (Ω). We recall that the space BV (Ω) of bounded variation functions (see [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF][START_REF] Aubert | Mathematical Problems in Image Processing, Partial Differential Equations and the Calculus of Variations[END_REF][START_REF] Attouch | Variational analysis in Sobolev and BV spaces : applications to PDEs and optimization[END_REF]) is defined as

BV (Ω) = {u ∈ L 1 (Ω) | Φ(u) < +∞},
where

Φ(u) = sup {∫ Ω u(x) div ξ(x) dx | ξ ∈ C 1 c (Ω), ξ ∞ 1 } . ( 1
)
The space BV (Ω), endowed with the norm u BV (Ω) = u L 1 + Φ(u), is a Banach space. The derivative in the sense of the distributions of every u ∈ BV (Ω) is a bounded Radon measure, denoted Du, and Φ(u) = ∫ Ω |Du| is the total variation of Du. We extend this definition to the second derivative (in the distributional sense). Recall that the Sobolev space

W 1,1 (Ω) = { u ∈ L 1 (Ω) | ∇u ∈ L 1 (Ω) }
where ∇u stands for the first order derivative of u (in the sense of distributions).

Definition 1. A function u ∈ W 1,1 (Ω) is Hessian bounded if |u| BV 2 (Ω) := sup {∫ Ω ∇u, div(φ) R n | φ ∈ C 2 c (Ω, R n×n ), φ ∞ 1 } < ∞, where div(φ) = (div(φ 1 ), div(φ 2 ), . . . , div(φ n )), with ∀i, φ i = (φ 1 i , φ 2 i , . . . , φ n i ) ∈ R n and div(φ i ) = n ∑ k=1 ∂φ k i ∂x k .
For more information on the BV 2 (Ω) space, see [START_REF] Demengel | Fonctions à hessien borné[END_REF][START_REF] Piffet | Modèles variationnels pour l'extraction de textures 2D[END_REF].

The variational model

We now assume that the image we want to recover from the data u d can be written as f = u + v where . where u is in BV (Ω) and v is in BV 2 (Ω). Such decompositions have already been performed [START_REF] Aubert | Modeling very oscillating signals, application to image processing[END_REF][START_REF] Aubert | Image decomposition into a bounded variation component and an oscillating component[END_REF][START_REF] Aubert | Mathematical Problems in Image Processing, Partial Differential Equations and the Calculus of Variations[END_REF] using the "Meyer" space of oscillating function [START_REF] Meyer | Oscillating patterns in image processing and nonlinear evolution equations[END_REF] instead of BV 2 (Ω). So far, the modelization we propose is not the same: the oscillating component will be included in the non regular part u while v involves the cartoon and the contours. We consider the following function defined on BV (Ω) × BV 2 (Ω) :

F (u, v) = 1 2 u d -u -v 2 L 2 (Ω) + λ|u| BV (Ω) + µ|v| BV 2 (Ω) + δ ∇v W 1,1 (Ω) , (2)
where λ, µ, δ ≥ 0. We are looking for a solution to the optimisation problem inf

(u,v)∈BV (Ω)×BV 2 (Ω) F (u, v) (P)
The first term

u d -u -v 2 L 2 (Ω)
of F is the fitting data term. Other terms are Tychonov-like regularization terms. Note that the δ-term is not useful from the modelling point of view. It is only a tool that allows to prove existence of solutions. We shall choose δ = 0 for numerical tests.

If the image is noisy, the noise is considered as a texture and will be involved in u: more precisely v will be the part of the image without the oscillating component, that is the denoised part. In [START_REF] Bergounioux | A BV 2 (Ω) model for image denoising and/or texture extraction[END_REF] we focus on the denoising process taking only v into account (and assuming that u = 0 so that u d -v is the noise). We now give an existence and uniqueness result for the general problem (P) (see [START_REF] Bergounioux | A BV 2 (Ω) model for image denoising and/or texture extraction[END_REF] for the proof).

Theorem 1. Assume that λ > 0, µ > 0 and δ > 0. Problem (P) has a unique solution (u, v).

Numerical implementation

Discretization of the problem

We assume for simplicity that the image is squared with size N × N . We note X := R N ×N R N 2 endowed with the usual inner product and the associated euclidean norm

u, v X := ∑ 1 i,j N u i,j v i,j , u X := √ ∑ 1 i,j N u 2 i,j . (3)
It is classical to define the discrete total variation as following (see for example [START_REF] Aubert | Mathematical Problems in Image Processing, Partial Differential Equations and the Calculus of Variations[END_REF]) : the discrete gradient of the numerical image u ∈ X is ∇u ∈ X 2 defined by (∇u) i,j = ( (∇u)

1 i,j , (∇u) 2 i,j ) , (4)
where (∇u)

1 i,j = { u i+1,j -u i,j if i < N 0 if i = N, and (∇u)
2 i,j = { u i,j+1 -u i,j if j < N 0 if j = N.
The (discrete) total variation |u| BV (Ω) is given by

J 1 (u) = ∑ 1 i,j N (∇u) i,j R 2 , (5)
where

(∇u) i,j R 2 = ( (∇u) 1 i,j , (∇u) 2 i,j ) R 2 = √ ( (∇u) 1 i,j ) 2 + ( (∇u) 2 i,j ) 2 .
The discrete divergence operator div is the adjoint operator of the gradient operator ∇ :

∀(p, u) ∈ X 2 × X, -div p, u X = p, ∇u Y , so that (div p) i,j =    p 1 i,j -p 1 i-1,j if 1 < i < N p 1 i,j if i = 1 -p 1 i-1,j if i = N +    p 1 i,j -p 2 i,j-1 if 1 < j < N p 2 i,j if j = 1 -p 1 i,j-1 if i = N. (6)
To define a discrete version of the second order total variation we have to introduce the discrete Hessian operator. As for the gadient operator, we define it by finite deferences. So, for any v ∈ X, the Hessian matrix of v, denoted Hv is identified to a X 4 vector: (Hv) i,j = ( (Hv) 11 i,j , (Hv) 12 i,j , (Hv) 21 i,j , (Hv) 22 i,j

) .

The discrete second order total variation |v| BV 2 (Ω) of v is defined as

J 2 (v) = ∑ 1 i,j N (Hv) i,j R 4 . (7)
As in the BV case, we may compute the adjoint operator of H (which is the dicretized "second divergence" operator) :

∀p ∈ X 4 , ∀v ∈ X H * p, v X = p, Hv X 4 . ( 8
)
and we deduce an numerical expression for H * from the equality [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. The discretized problem stands inf

(u,v)∈X 2 1 2 u d -u -v 2 X + λJ 1 (u) + µJ 2 (v) + δ(|v| + J 1 (v)), (P d )
where

|v| := ∑ 1 i,j N |v i,j | .
In the finite dimensional case we have an existence result with δ = 0.

Theorem 2. Problem P d has a unique solution for every δ ≥ 0 and λ > 0.

For numerical purpose we shall set δ = 0. In fact, we have performed tests with δ = 0 and very small δ = 0 (as required by the theory to get a solution to problem P) and results where identical. So, to simplify numerical implementation, we consider the following discretized problem :

inf (u,v)∈X 2 1 2 u d -u -v 2 X + λJ 1 (u) + µJ 2 (v). ( Pd )

Algorithm

Using non smooth analysis tools (for convex functions) it is easy to derive (necessary and sufficient) optimality conditions. More precisely (u, v) is solution of ( Pd ) if and only if

   0 ∈ ∂ ( λJ 1 (u) + 1 2 u d -u -v 2 ) 0 ∈ ∂ ( µJ 2 (v) + 1 2 u d -u -v 2 ) , ( 9
)
where ∂J is the classical subdifferential of J. Using subdifferential properties, we see that ( 9) is equivalent to

{ u = u d -v -Π λK1 (u d -v) v = u d -u -Π µK2 (u d -u) . ( 10
)
where K 1 and K 2 are closed convex sets. Chambolle [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] proved that

K 1 = {div p | p ∈ Y, p i,j R 2 1 ∀i, j = 1, . . . , N } (11)
in the BV (Ω) setting and we may prove similarly that

K 2 = {H * p | p ∈ Z, p i,j R 4 1, ∀i, j = 1, . . . , N }, ( 12
)
(see [START_REF] Bergounioux | A BV 2 (Ω) model for image denoising and/or texture extraction[END_REF]). Moreover, Chambolle [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] proposed a fixed point algorithm to compute Π λK1 and we are able to extend this result to the second order case. Therefore, we set the following fixed point algorithm:

p 0 = 0 (13a) p n+1 i,j = p n i,j -τ (H[H * p n -u d /λ]) i,j 1 + τ (H[H * p n -u d /λ]) i,j R 4 . ( 13b
)
which convergence is proved in [START_REF] Bergounioux | A BV 2 (Ω) model for image denoising and/or texture extraction[END_REF] :

Theorem 3. Let be τ 1/64. Then λ (H * p n ) n converges to P λK2 (u d ).
So, we propose the following algorithm :

• Step 1 : We choose u 0 et v 0 (for example, u 0 = 0 et v 0 = u d ) and 0 < α < 1.

• Step 2 : define the sequences ((u n , v n )) n as follows:

{ u n+1 = u n + α (u d -v n -Π λK1 (u d -v n ) -u n ) v n+1 = v n + α (u d -u n -Π µK2 (u d -u n ) -v n ) .
• Step 3 : if a stopping criterion is not satisfied, set k := k + 1 and go back to 2.

We are not able to prove the convergence of this algorithm, but we observed that, numerically, it always converges for appropriate values of α. Nevertheless, we propose in [START_REF] Piffet | Modèles variationnels pour l'extraction de textures 2D[END_REF] another algorithm which convergence is rigorously proved .

Numerical tests and comments

We test the model on two images: The first one is a synthetic image wher texture has been artificially added, and the second one is the well known "Barbara" benchmark, often used in texture extraction. We perform many tests with respect to the different parameters. We only present here the most significant : α = 0.6, λ = 0.5 and µ = 100. Let us first report on the iterations number effect with image (a).

If we are only interested in the texture part, we can observe in fig 2 that we get back all the texture. Unfortunately, most of geometrical information (that we don't want to get back here) is also kept, and we observe the involved geometric part is getting more important as the iteration number is growing. In addition, we see that the too many geometrical information remains together with the texture in the oscillating part: this is a bad point. Nevertheless, our main goal is to locate the texture and we don't need to work with the cartoon part anymore once it has been identified. We do not need to recover all the texture but only a significant part to identify it. In that case, we propose a method that permits to improve the results significantly: we modify the hessian operator to make it anisotropic. More precisely, we reinforce chosen directions. As texture is made of oscillating information, we hope that we shall keep most of it while many contour lines disappear. We specially, act on the vertical and horizontal components of the hessian operator. To deal with non vertical and horizontal lines, we just have to let the image rotate. In the following test, we have replaced the hessian operator by the operator H defined for all v ∈ X 4 by : ∀(i, j) ∈ {1, ..., N } 2 , (H v) i,j = ( 0, (Hv) 12 i,j , (Hv) 21 i,j , 0

) .

Conclusion

The model permits to extract texture from an image, but the texture part still contains too much geometric information. Thus, to recover what we are interested in, we have to use the algorithm with a limited number of iterations. Moreover, we have noticed that we recover too many contour lines as well. The asset of this model is that we can make it anisotropic, modifying the hessian operator in an appropriate way. Therefore we get rid of geometrical information, but we lose part of the texture as well. Nevertheless, if our goal is just to locate texture on an image, this loss remains acceptable.

Fig. 1 .

 1 Fig. 1. Original images.

Fig. 2 .

 2 Fig. 2. Number of iterations: first line: 60; second line: 200; third line: 2000.

Fig. 3 .

 3 Fig. 3. Number of iterations: 2000.

Fig. 4 .

 4 Fig. 4. Test with the anisotropic operator H . Number of iterations: first line: 60; second line: 2000.


Bart Lamiroy

Jean-Philippe Ropers

Assessing Inductive Logic Programming Classification Quality by Image Synthesis

One of the major difficulties in classification is the assessment of classification quality. This paper builds on a recently developed classification technique, based on Inductive Logic Programming to propose a way of visually reconstructing the learnt classes in order to give a feedback to the user. The system translates positioning constraints into a linear programming problem, which can be solved with standard state-of-theart approaches. The result can then be used for generating a synthetic image, representative of the class.

Introduction

We have recently [1,[START_REF] Lamiroy | Utilisation de Programmation Logique Inductive pour la reconnaissance de symboles[END_REF] introduced a new approach to graphical symbol representation and classification by means of a restricted visual vocabulary and the use of Inductive Logic Programming (ILP) [3]. Our current investigations have shown that ILP is a very powerful framework and that it allows for types of symbol characterization of symbols that are not possible with numerical classification techniques. Since one of the advantages of ILP is to be able to express image properties in a close-to natural language, it is necessary to have an easy, but objective way to assess the image properties that are used for classification.

In this paper we develop a method that allows us to generate a synthetic image from a set of first order logic expressions that are characterizing an image class. This will allow the user to validate if the learning process has been coherent with the context semantics of the images under consideration.

Visual Vocabulary and ILP

The approach developed in [1,[START_REF] Lamiroy | Utilisation de Programmation Logique Inductive pour la reconnaissance de symboles[END_REF] and recently extended, produces a description of a symbol using first order logic predicates. These predicates express relations between previously extracted visual cues from the image (composing the visual vocabulary). These cues are extracted using classical image treatment algorithms and consist of: circles [START_REF] Lamiroy | Robust Circle Detection[END_REF], loose endpoints, thick components and corners.

These cues are then linked together with relative positioning relations (north, north-west, west ...), bigger, smaller, near, far, inside, superimposed. As a result, the image description can be expressed with first order logic predicates that are quite close to a natural language description like: "Image A contains two circles, one upper-right corner and two lower-left corners; the first circle lies left (east) to the second, and above (north) the upper-right corner ...". More prosaically speaking, this gives something like: type(prim_A,circle). type(prim_B,cornerne). type(prim_C,cornersw). type(prim_D,cornersw). type(prim_E,circle). has_element(img,prim_A). has_element(img,prim_B). has_element(img,prim_C). has_element(img,prim_D). w(prim_A,prim_E). n(prim_B,prim_A). ... Inductive Logic Programming allows to apply supervised learning on these data, provided that it's given 1. a set of known vocabulary, rules, axioms or predicates, describing the domain knowledge base K; 2. a set of positive examples E + the system is supposed to describe or characterize with the set of predicates of K ;

3. a set of negative examples E -that should be excluded from the deducted description or characterization.

The ILP solver is then able to find the set of properties P, expressed with the predicates and terminal vocabulary of K such that the largest possible subset of E + verifies P, and such that the largest possible subset of E -does not verify P. This approach has already been successfully used for extraction of semantics from written text [START_REF] Claveau | From efficiency to portability: Acquisition of semantic relations by semisupervised machine learning[END_REF] or in document and image analysis structures [START_REF] Amin | Learning to recognize hand-printed chinese characters using inductive logic programming[END_REF][START_REF] Ceci | Relational data mining and ilp for document image processing[END_REF]. This approach allows for learning common properties within classes of symbols such as to express non-trivial knowledge of visual representation of more semantic concepts.

The conclusions of the previous studies have shown that ILP is very useful for describing complex and non trivial symbols and taking into account non-geometric transformations. However, the approach remains tributary of the initially defined vocabulary K. More specifically the adjunction of new relationships influences on the classification "quality" (provided there is such a notion of quality). In this paper we propose a subjective classification assessment tool, based on image synthesis.

The paper is organized as follows. First we give a brief overview of the classification differences induces by variations in the description vocabulary (section 2.2). We then provide an approach to retrieve image geometry from an first order logic image description in section 3, before concluding and opening new perspectives in section 2.4.

Vocabulary Description

In order to give a clear example of what can vary in the classification results, depending on the used vocabulary, we introduce different types of predicates. Relations like proximity, size etc.. significantly modify classification results. In this section we describe the semantics of the new relations we are adding to K. These new relations concern: connectedness, qualitative distance and relative size [START_REF] Mas Romeu | An Incremental On-line Parsing Algorithm for Recognizing Sketching Diagrams[END_REF]. We shall then combine these in order to see how they influence on the final classification.

Inclusion and Connectedness

In order to distinguish whether two elements in the image are connected we construct the relation superimposed and set up a new algorithm to evaluate if two elements are bound. To do that, we compute the connected component tree of each individual element and the connected component tree of synthetic image summing both. Since, by construction of the different vocabulary-related extraction algorithms, each element has a unique black component at the root of this tree, it is quite straightforward to determine if two components are connected, disjoint or included one into another by simply counting the number of components of the fusion of both elements. Indeed, if the number of connected components of the fusion is equal to the sum of connected components of each individual element, then both are disconnected. On the contrary, if the depth of the resulting component tree is equal to the sum of individual depths, then both are included. Any other situation means both elements are connected (and therefore superimposed).

Dimensional Relations

The implementation of the dimensional relations are also based on the connected component trees. We proceed by comparing the area of the first child of the root of both elements. All elements are linked pairwise by a size relation : bigger, smaller or samesize. Here again, we used the hypothesis that the connected component tree of a descriptor is a root with only one child.

This algorithm is based on a pixel-wise estimation of the area covered by the element, and is therefore a very reliable measure.

Proximity

While the two previous sets of relations were quite obvious, proximity is dealt with in a more subtle manner. near and far are subjective (or at least non-absolute) concepts. On our case, these relations are considered to be "object-centric", meaning that B is close with respect to A if the distance between both is comparable to the size of A. This induces, however, that the proximity relation is asymmetric.

As a practical point of view, we compute first a dilation of the element A, then if A and B are superimposed (according to relation above) B is close to A. The size of the structural element is

Influences on Classification

In what follows we are going to compare different configurations and see how changes in the vocabulary affect learned classification rules. We distinguish five configurations:

• the basic configuration, including directional relations with inside as in [START_REF] Lamiroy | Utilisation de Programmation Logique Inductive pour la reconnaissance de symboles[END_REF];

• the proximity configuration, directional relations, near and far as described in section 2.1.3;

• the superimpose configuration: directional relations, superimposed and inside described in section 2.1.1;

• the size configuration: directional relations, bigger, smaller, samesize as described in section 2.1.2;

• all the relations: directional relations, near, far, superimposed, inside, bigger, smaller and samesize.

Introducing Ambiguity

We use the set of symbols from Figure 2 and use our method to learn a description of these symbols with respect to a set of counter-examples. The different configurations (i.e. allowing more or less predicates to be used, as mentioned above) we have tested give the following, results:

• basic: has_element(A,B), type(B,circle), has_element(A,C), inside(C,B).

• proximity: has_element(A,B), type(B,circle), has_element(A,C), near(B,C).

• superimpose: has_element(A,B), type(B,circle), has_element(A,C), superimposed(C,B).

• size: has_element(A,B), type(B,circle), has_element(A,C), bigger(B,C).

• all: has_element(A,B), type(B,circle), has_element(A,C), bigger(B,C).

Although this might seem normal at a first glance (actually it is quite comforting as a result) it raises a fundamental problem: the one of the uniqueness of the symbol representation. All configurations give perfect classifications, and use the same number of predicates and of primitives. Although one can argue that the first four situations cannot be compared per se, the last one creates the issue. What causes our system to select bigger as distinguishable predicate, compared to superimposed, near or inside? This needs to be further investigated, and necessarily depends on the implementation, order of declaration or any other arbitrary rule. It is clear that there will be a strong need to define an order of priority on our predicate ruleset, at least. Although this is essentially not a problem with, for instance, SVM classifiers [9], it raises an interesting question on classification. Indeed, all classification approaches are evaluated on their capacity of correctly separating instances from one another and affect a class label to them and are tested with respect to a closed world of instances. The fact that they are usually trained on a distinct training set, and tested on another set doesn't quite matter that much: there is no real evaluation on the sensitivity of the classification boundaries with respect to the training set, as long as the final classification on the test set is acceptable. The fact that our approach, compared to others (i.e. statistical) gives descriptions that are readily interpretable allows for a human-based validation of the classification rules, and can thus be embedded in a relevance feedback process.

Shorter Rules

In most experiments, the fact of adding new predicates results in shorter classification clauses, as shown on the example below. The basic configuration, expresses that the diode symbols represented in Figure 3 With the addition of the newly defined relations, the rule produced is one element shorter.

symbol(A) :-has_element(A,B), type(B,ccblackthick), has_element(A,C), has_element(A,D), e(D,C), superimposed(D,B).

We can see that superimposed replaces inside and a cornerse.

Adding Complexity

Our final series of experiments concern situations where the positive example configurations cannot be reduced to a single predicate. The fact of adding predicates does not fundamentally increase the number classification rules. The rules produced this time are a bit different but the structure is almost the same. One of the rule sets produced in the basic case is: Here again, as in the previous section, it would be an error to assess this result from a human-centric point of view. The fact that two rules describing the symbol class with respect to the presence of cornerne and cornerse elements suddenly is described by the presence of a cornersw element is only due to the fact that our approach is just another minimization technique, ILP looks for the most generic ruleset. This means containing the smallest number of predicates, covering the largest set of positive examples and avoiding the largest set of negative examples.

[Rule 1] [Pos cover = 1 Neg cover = 0] symbol(img_180_1).

Conclusion on Classification

The conclusion of this work is in line with previous work [START_REF] Lamiroy | Utilisation de Programmation Logique Inductive pour la reconnaissance de symboles[END_REF] and with what we expected: the enhanced expressiveness of the description language improves the produced classification rules, by reducing their size. It would be an error to consider that the increase of the vocabulary would induce a richer description of the classified symbols. As mentioned in section 2.3, this would be contrary to the optimization criteria of our approach, which tries to minimize the number of predicates used to describe a symbol.

One important issue has not been covered by this work, is the study of the the effect of the vocabulary increase with respect to the method's inherent sensitivity to the detection of the visual elements and their associated uncertainty. One of the main bottlenecks is the fact that ILP doesn't allow uncertainty to be integrated in its ruleset. One of the extensions of the method presented in this paper will be the use of an uncertainty-aware solver [START_REF] Richardson | Markov logic networks[END_REF][START_REF] Domingos | Markov logic[END_REF].

On the other hand, it is clear that the final results may seem awkward with respect to a human-centric point of view. The rest of this paper will therefore develop a method allowing to represent the obtained first-order logic description into a visual form, in order to visually asses the sense of the obtained set of predicates.

Generating Images

Once we have a description of an image in our previously described language, it is quite straightforward to recover an image that is coherent with this description. In order to achieve this, we simply express the relations as a Linear Programming [START_REF] Dantzig | Linear Programming and Extensions[END_REF] problem.

Modeling the Problem

As a matter of fact, the set of predicates describing an image can be seen as a set of constraints over the properties of all components composing the image, and described in the previous sections. As we shall show in section 3.2, these constraints can all be approximated by linear inequations. Combining this set of linear inequations with an objective function f expressing that the resulting arrangement of components should be as compact as possible, we obtain a linear constraint set A and a linear objective function f such that, with

x = (x 1 , x 2 , . . . x n ), c = (c 1 , c 2 , . . . c n ) and f (x 1 , x 2 , . . . x n ) = c 1 x 1 +c 2 x 2 +. . . c n x n , we obtain f (x) = xc t
and Ax ≤ b for minimization.

Inequation Set

In order to represent our image description with linear constraints, we are going to model all components by their bounding box. Let C be a detected component. The centre of its bounding box is given by (x C , y C ), its width and height by w C and h C .

Represented in the previously described linear programming framework, this gives for n visual components included in an image I.

x = (x C 1 , y C 1 , h C 1 ,
For each identified type of component, extra constraints apply: for circles, w C = h C ; extremities have w C = h C = 1 and corners have w C = h C = t where t is the template size used for detection.

The relationships between any two components C 1 and C 2 give the following constraints. It is to be noted that, due to our deliberate choice of using linear programming, the constraints are reduced to their linear expression:

• bigger(C 1 ,C 2 ) gives w C 1 -w C 2 ≥ 0 and h C 1 -h C 2 ≥ 0.
smaller gives rise to the dual equations.

• inside(C 1 ,C 2 ), needs a little quirk to get linearised. Indeed, the condition for bounding box C 1 being entirely included into C 2 is the following:

w C 2 > w C 1 (1) |x C 1 -x C 2 | + w C 1 2 < w C 2 2 (2)
(and similarly for y and h parameters). Unfortunately, the absolute value function is non linear. However, we can constrain our images to have positive coordinates. Under this assumption, the above constraints can be reformulated as

w C 2 > w C 1 (3) x C 1 -x C 2 + w C 1 2 < w C 2 2 (4) x C 2 -x C 1 + w C 1 2 < w C 2 2 (5)
which are perfectly linear.

• superimposed(C 1 ,C 2 ) is approximated by the overlapping of bounding boxes:

(x C 1 , w C 1 ) and (x C 2 , w C 2 ) overlap iff |x C 1 -x C 2 | ≤ w C 1 +w C 2 2
. We can however get away with the same linearity trick as before such that superimposed(C 1 ,C 2 ) gives

x C 1 -x C 2 ≤ w C 1 + w C 2 2 (6) x C 2 -x C 1 ≤ w C 1 + w C 2 2 (7)
(and similarly for y and h parameters).

• near(C 1 ,C 2 ) is more subtle to model linearly, since it is based on a non-linear distance metric. As described in section 2.1.3, proximity is relative to the size of the first argument

C 1 : C 2 is near to C 1 , if it is overlaps with a dilation of C 1 with a structural element of size w C 1 × h C 1 .
As a consequence, near can be expressed using the same constraints as the superiposed relation, by using

w C 1 = 2 × w C 1 and h C 1 = 2 × h C 1 .
• far(C 1 ,C 2 ) is, unfortunately not appropriately representable in our framework. There is no straightforward, linear way to express the fact that two elements are either not near, or explicitly relate the parameters as to model their distance with respect to a given threshold. The linearisation trick used for the near, superimposed or inside predicates cannot be applied here, for the reason that we need to express a ≥ relationship, rather than a ≤ as it was the case before.

The solution comes from the objective function we will discuss in the following paragraph.

The inequation set obtained from the set of predicates, combined with the objective function to minimize is then simply given to a linear programming constraint solver [START_REF] Dantzig | Linear Programming and Extensions[END_REF] making it very straightforward to generate the corresponding image. In order to express as closely as possible the fact that non-near elements should be as far away from each other as possible, we use an objective function that maximizes the dispersion of the elements. Here again, dispersion is a non-linear measure, which we need to approximate.

In our case, the adopted solution is to maximize

f (x) =   ∑ C i =I (x C i -x I ) + (y C i -y I )   -   ∑ C i =I w C i + h C i   -w I -h I
This will tend to minimize the size of the components while "maximizing' their "distance' from the centre of the image.

Experimental Results

In order to assess the validity of our approach we first try to synthesize the classification results obtained in section 2.2. The classification rules are given below, followed by the generated, corresponding image.

Figure 1 :

 1 Figure 1: relations far and near

Figure 2 :

 2 Figure 2: Image Set of Positive Examples.

  contain a ccblackthick and a cornerse component, as well as an unqualified component D. The corner is east of D and the thick component is inside D.

Figure 3 :

 3 Figure 3: Diode Examples.

[

  Rule 2] [Pos cover = 2 Neg cover = 0] symbol(A) :-has_element(A,B), type(B,circle), has_element(A,C), e(B,C), type(C,cornerne). [Rule 3] [Pos cover = 1 Neg cover = 0] symbol(img_184_1). [Rule 4] [Pos cover = 2 Neg cover = 0] symbol(A) :-has_element(A,B), type(B,ccblackthick), has_element(A,C), has_element(A,D), e(B,C), type(C,cornerse), inside(D,B). Introducing the new predicates changes the ruleset to three individual rules and one collective covering three examples, this rule is depicted below: symbol(A) :-has_element(A,B), type(B,ccblackthick), has_element(A,C), has_element(A,D), e(C,D), far(D,C), type(D,cornersw).

  w C 1 , . . . x Cn , y Cn , h Cn , w Cn , x I , y I , h I , w I )

ConclusionsIn this paper, the components in an image are represented by their bounding box. This might induce errors with complex shapes. Without loss of generality, components can be approximated by more complex boundaries, while remaining linear. It will have to be established what impact this might have on time complexity, however. Secondly, the relations between components are approximated with linear equations, while they are essentially non-linear. This approximation is limiting for more complex shape descriptions. Further work needs to be done to exactly evaluate the impact of this approximation in terms of speed, found optimum or possible singularities and the authors are currently working on a non-linear extension of the minimization algorithm.

Acknowledgements

Part of the experiments described in this paper are due to the Master Project of Karly Langa and Benoît Leoutre, students at ESIAL -Nancy-Université.