Abdurrahim Akgundogdu
email: akgundog@istanbul.edu.tr

Rachid Jennane

Gabriel Aufort

Claude Laurent Benhamou

Osman Nuri Ucan

3D IMAGE ANALYSIS AND ARTIFICIAL INTELLIGENCE FOR BONE DISEASE CLASSIFICATION

Keywords: Trabecular bone, Hybrid Skeleton Graph Analysis, SVM, GA, Anfis

In order to prevent bone fractures due to disease and ageing of the population, and to detect problems while still in their early stages, 3D bone micro architecture needs to be investigated and characterized.

Here, we have developed various image processing and simulation techniques to investigate bone micro architecture and its mechanical stiffness. We have evaluated morphological, topological and mechanical bone features using artificial intelligence methods. A clinical study is carried out on 2 populations of arthritic and osteoporotic bone samples. The performances of Adaptive Neuro Fuzzy Inference System (ANFIS), Support Vector Machines (SVM) and Genetic Algorithm (GA) in classifying the different samples have been compared. Results show that the best separation success (100 %) is achieved with Genetic Algorithm.

Introduction

The interest in three-dimensional image analysis is increasing in different fields of image processing. In biomedical systems many new tomographic modalities now model various 3-D images of organs and bones. Applications of these methods range from biological studies [START_REF] Thiel | Shape splitting from medial lines using the 3-4 chamfer distance[END_REF] to character recognition [START_REF] Gonzalez | Digital Image Processing[END_REF]. Image analysis techniques are used in order to obtain morphological and topological characterizations.

In the case of bone, the micro architecture plays an important role in osteoporosis, osteoarthritis and the prediction of fractures. To characterize such material, Cruz-Orive [START_REF] Cruz-Orive | Estimation of sheet thickness distrib ution from linear and plate Sections[END_REF] and Levitz et al. [START_REF] Levitz | Disordered porous solids: from chord distributions to small angle scattering[END_REF] evaluated the mean size of materials by averaging 3D data. Odgaard et al. [START_REF] Odgaard | Quantification of connectivity in cancellous bone, with special emphasis on 3D reconstructions[END_REF] and Vogel [START_REF] Vogel | Digital unbiased estimation of the Euler-Poincaré characteristic in different Dimensions[END_REF] defined a connectivity index that can be computed by scanning data and storing topological information. In addition to these methods, the morphological structure can be accessed by determining parameters that can be computed directly from the 3D image without an underlying model assumption [START_REF] Hildebrand | A new method for the model-independent assessment of thickness in three-dimensional images[END_REF][START_REF] Jennane | Fractal Analysis of Bone X-Ray Tomographic Microscopy Projections[END_REF]. This work led however to global methods based on physical models that cannot give precise information about the medium's structure and its local properties. Nowadays, it is possible to study the object at a local scale and obtain information from each of the elements that compose the structure of the material [START_REF] Pothuaud | A new method for three-dimensional skeleton graph analysis of porous media : application to trabecular bone microarchitecture[END_REF]. A skeleton-based technique called Line Skeleton Graph Analysis (LSGA) [START_REF] Aufort | A New 3D Shape-Dependent Skeletonization Method[END_REF] enables the bone micro architecture to be studied at a local scale, providing global information about the bone. This method preserves the topology of the medium but is not suitable for structures comprising different shapes such as trabecular bone, which is composed of rod-shaped and plate-shaped elements. In this case, all non-cylindrical shapes are better described by 2D-surfaces rather than by 1D-curves.

We have developed a new technique called Hybrid Skeleton Graph Analysis (HSGA) [START_REF] Aufort | Hybrid skeleton graph analysis of disordered porous media. Application to trabecular bone[END_REF] to create structural models of disordered porous media. The method is based on [START_REF] Pothuaud | A new method for three-dimensional skeleton graph analysis of porous media : application to trabecular bone microarchitecture[END_REF], but takes into account the shape of the object, in order to improve the geometrical approximation of [START_REF] Pothuaud | A new method for three-dimensional skeleton graph analysis of porous media : application to trabecular bone microarchitecture[END_REF], since disordered porous media such as trabecular bone are not homogeneous but are composed of rodshaped and plate-shaped elements. Curve thinning and surface thinning are efficient respectively for one or the other, but not for both. Our early work for skeletonizing hybrid-shaped media focused on providing a tool for objects composed of rods and plates [START_REF] Aufort | A New 3D Shape-Dependent Skeletonization Method[END_REF]. The HSGA method is based on a homotopic hybrid skeleton using curve and surface thinning techniques which preserve Betti-numbers [START_REF] Morgenthaler | Three-dimensional simple points: serial erosion, parallel thinning and skeletonization[END_REF]; the topology of the medium is thus perfectly preserved.

At the same time, the combination of skeletons and Finite Elements (FE) analysis to evaluate the stiffness of large-scale porous media has recently shown great potential [START_REF] Pothuaud | A new computational efficient approach for trabecular bone analysis using beam models generated with skeletonized graph technique[END_REF][START_REF] Lenthe | Specimen specific beam models for fast and accurate prediction of human trabecular bone mechanica l properties[END_REF]. We have improved on a previously published beam Finite Element Analysis (FEA) in [START_REF] Pothuaud | A new computational efficient approach for trabecular bone analysis using beam models generated with skeletonized graph technique[END_REF]. A full protocol [START_REF] Aufort | Nouvelle approche de modélisation de milieux poreux. Application à l'os trabéculaire[END_REF][START_REF] Aufort | Mechanical assessment of porous media using hybrid skeleton graph analysis and finite elements. Application to trabecular bone[END_REF] has been developed for evaluating the stiffness of large scale porous media also based on the HSGA and using the advantage of both beam and shell Finite Elements to enable a fast and precise mechanical simulation.

The improvement of these techniques has enabled the determination of various parameters for quantifying the morphology and the stiffness of the bone micro architecture.

These parameters have to be combined efficiently to study bone diseases. For this purpose, knowledge modeling is achieved using Artificial Intelligence (AI) approaches which are developed by human expertise [START_REF] Cherkassky | Fuzzy Inference Systems: A Critical Review, Computational Intelligence, Soft Computing and Fuzzy-Neuro Integration with Applications of NATO ASI Series[END_REF][START_REF] Jang | Self-learning fuzzy controllers based on temporal backpropagation[END_REF][START_REF] Vapnik | The nature of statistical learning theory[END_REF][START_REF] Siedlecki | A Note on Genetic Algorithms for Large-Scale Feature Selection[END_REF]. These models are based on training and testing. The Adaptive Neuro Fuzzy Inference System (ANFIS), Support Vector Machines (SVM) and Genetic Algorithm (GA) are the most successful tools for predicting input-output data set systems.

This work presents a study performed on 2 populations of trabecular bone samples with an a priori knowledge of fracture risk, 9 osteoarthritic (OA) samples and 9 osteoporotic (OP) samples. The different improved techniques have been processed on each of the 18 bone samples, producing features which are then used to quantify the microarchitectural bone quality. These parameters are used as entries for different artificial intelligence systems, namely: Adaptive Neuro-Fuzzy Inference System (ANFIS), Support Vector Machines (SVM) and Genetic Algorithm (GA). We also investigated whether the combination of these different techniques contributes to a better prediction of the nature (OA or OP) o f each sample studied.

The paper is organized as follows: first, the different artificial intelligence systems used (ANFIS, SVM and GA) are presented. Then, the techniques used to compute the OA and OP data are presented. The results from the clinical study follow. A discussion concludes this work.

Material and Methods

This section presents the three artificial intelligence methods, Adaptive Neuro Fuzzy Inference System (ANFIS), Support Vector Machines (SVM) and Genetic Algorithm (GA), as well as the image processing techniques, Hybrid Skeleton Graph Analysis (HSGA) and Finite Element Analysis (FEA) used to classify each sample of the two studied populations (OA and OP).

Artificial Analysis

Artificial Intelligence (AI) is the area of computer science which deals with machine help in finding solutions to complex problems in a more human-like fashion. This generally involves borrowing characteristics from human intelligence, and applying them as algorithms in a computer-friendly way. ANFIS is a combination of ANN and Fuzzy Inference System (FIS) [START_REF] Jang | Self-learning fuzzy controllers based on temporal backpropagation[END_REF] which is used to determine the parameters of FIS. ANFIS implements a Takagi Sugeno Kang (TSK) [START_REF] Sugeno | Structure identification of fuzzy model[END_REF] fuzzy inference system in which the conclusion of the fuzzy rule is constituted by a weighted linear combination of the crisp inputs.

Architecture of ANFIS

ANFIS is a fuzzy Sugeno model implemented in the framework of adaptive systems to facilitate learning and adaptation [START_REF] Cherkassky | Fuzzy Inference Systems: A Critical Review, Computational Intelligence, Soft Computing and Fuzzy-Neuro Integration with Applications of NATO ASI Series[END_REF]. Such a framework makes ANFIS modeling more systematic and less reliant on expert knowledge. To present the ANFIS architecture, two fuzzy if-then rules based on a first order Sugeno model are considered [START_REF] Cherkassky | Fuzzy Inference Systems: A Critical Review, Computational Intelligence, Soft Computing and Fuzzy-Neuro Integration with Applications of NATO ASI Series[END_REF][START_REF] Jang | Self-learning fuzzy controllers based on temporal backpropagation[END_REF][START_REF] Sugeno | Structure identification of fuzzy model[END_REF][START_REF] Sun | Rulebase structure identification in an adaptive-network-based fuzzy inference system[END_REF].

Rule 1: If (x is A 1 ) and (y is B 1 ) then 1 1 1 1 f p q y r    , (1)
Rule 2: If (x is A 2 ) and (y is B 2 ) then 2 2 2 2 f p q y r    . (2)
where x and y are the inputs, A i and B i are the fuzzy sets, f i are the outputs within the fuzzy region specified by the fuzzy rule, and i i i p , q and r are the design parameters that are determined during the training process. The ANFIS architecture for the implementation of these two rules is shown in Figure 1, in which a circle indicates a fixed node, whereas a square indicates an adaptive node. In layer 1, all the nodes are adaptive nodes. The outputs of layer 1 are the fuzzy membership grade of the input, which are given by

1 i O = ( ) Ai x  , i= 1, 2, or (3)
1 i O i = 2 ( ) B i y   , i= 3, 4, (4)
where

() Ai x  , 2 () Bi y  
can adopt any fuzzy membership function. We will use the Gaussian membership function given by:  are parameters to be learnt. These are the premise parameters.

In layer 2, the nodes are fixed. They are labeled with M, indicating that they behave as a simple multiplier. The outputs of this layer can be represented as 2 ( ) ( ), 1, 2

i i i i A B O w x y i      (6)
which are the so-called firing strengths of the rules.

In layer 3, the nodes are also fixed nodes. They are labeled with N, indicating that they play a normalization role to the firing strengths from the previous layer. The outputs of this layer can be represented as

3 1 2 , 1, 2 i i i w O w i w w     (7)
These are the so-called normalized firing strengths.

In layer 4, the nodes are adaptive nodes. The output of each node in this layer is simply the product of the normalized firing strength and a first-order Sugeno model. Thus, the outputs of this layer are given by 4 ( ), 1, 2

i i i i i i i O w f w p x q y r i      (8)
In layer 5, there is only one single fixed node labeled with S. This node performs the summation of all incoming signals. Hence, the overall output of the model is given by

2 2 5 1 1 1 2 ( ) . i i i i i i i w f O w f w w        (9)
It can be observed that there are two adaptive layers in this ANFIS architecture, namely the first layer and the fourth layer. In the first layer, there are three modifiable parameters  

,, iii a b c , which are related to the input membership functions. These parameters are the so-called premise parameters.

In the fourth layer, there are also three modifiable parameters   ,,

i i i
pqr , pertaining to the first order Sugeno model. These parameters are called consequent parameters [START_REF] Sun | Rulebase structure identification in an adaptive-network-based fuzzy inference system[END_REF][START_REF] Takagi | Derivation of fuzzy control rules from human operator's control actions[END_REF].

Hybrid Learning Algorithm of ANFIS

The task of the learning algorithm for this architecture is to tune all the modifiable parameters, namely   ,,

iii a b c and   ,,

i i i
pqr, in order to make the ANFIS output match with the training data.

When the premise parameters ,,

iii a b c of the membership function are fixed, the output of the ANFIS model can be written as:

1 2 1 2 1 2 1 2 w w f f f w w w w     . ( 10
) 1 2 1 2 . f w f w f   (11) 1 2 1 1 1 2 2 2 ( ) ( ). f w p x q y r w p x q y r       ( 12
) 1 2 2 2 1 1 1 1 1 2 2 2 ( ) ( ) ( ) ( ) ( ) ( ) f w x p w y q w r w x p w y q w r       . ( 13
)
Where 1 1 1 2 2 2 , , , , , p q r p q r are linear consequent parameters. When the premise parameters are not fixed, the search space becomes larger and the convergence of the training becomes slower. A hybrid algorithm combining the least squares method and the gradient descent method is adopted to solve this problem [START_REF] Übeyl | Automatic detection of erythemato-squamous diseases using adaptive neuro-fuzzy inference systems[END_REF]. The hybrid algorithm is composed of forward pass and backward pass.

The least squares method (forward pass) is used to optimize the consequent parameters with the fixed premise parameters. [START_REF] Takagi | Fuzzy identification of systems and its applications to modeling and control[END_REF][START_REF] Fogelman | Neuro-computing: Algorithms, Architectures and Applications[END_REF]. Once the optimal consequent parameters have been found, the backward pass starts immediately. The gradient descent method (backward pass) is used to optimize the premise parameters corresponding to the fuzzy sets in the input domain. The output of the ANFIS is calculated by employing the consequent parameters found in the forward pass.

The output error is used to adapt the premise parameters by means of a standard back propagation algorithm. It has been proven that this hybrid algorithm is highly efficient in training the ANFIS [START_REF] Fogelman | Neuro-computing: Algorithms, Architectures and Applications[END_REF][START_REF] Broomhead | Multivariable functional interpolation and adaptive networks[END_REF].

Support Vector Machines (SVM)

Support Vector Machine is an effective technique for data classification. The SVM proposed by Vapnik [START_REF] Vapnik | The nature of statistical learning theory[END_REF] has been extensively studied for classification, regression and density estimation. A classification task usually involves training and testing data which consist of some data instances.

Each sample in the training set contains one "target value" (class labels) and several "attributes" (features). The goal of SVM is to produce a model which predicts the target value of data instances in the testing set where only the attributes are known. Figure 2 illustrates the map of inputoutput SVM space. In Figure 3, there are many possible linear classifiers that can separate the data, but there is only one that maximizes the margin (maximizes the distance between the nearest data point of each class). This linear classifier is termed the optimal separating hyper plane. Viewing the input data as two sets of vectors in an n-dimensional space, an SVM will construct a separating hyper plane in that space, which maximizes the "margin" between the two data sets. To calculate the margin, we construct two parallel hyper planes, one on each side of the separating one, which are "pushed up against" the two data sets. Intuitively, a good separation is achieved by the hyper plane that has the largest distance from the neighboring data points of both classes. It is expected that the larger the margin or distance between these parallel hyper planes, the better the generalization error of the classifier will be.

Let     , ,
k i i i i f x sign y K x x b              ( 15
) 1 0 sgn( ) 1 0 for u u for u       
where k is the number of learning patterns, i y is the target value of learning pattern The idea of the kernel function is to enable operations to be performed in the input space rather than the potentially high dimensional feature space. Hence the inner product does not need to be evaluated in the feature space. This provides a way of addressing the problem of dimensionality. However, the computation is still critically dependent upon the number of training patterns and providing a good data distribution for a high dimensional problem will generally require a large training set. This theory is based upon Reproducing Kernel Hilbert Spaces (RKHS) [START_REF] Aronszajn | Theory of reproducing kernels[END_REF][START_REF] Girosi | An equivalence between sparse approximation and Support Vector Machines[END_REF][START_REF] Wahba | Spline Models for Observational Data[END_REF]. An inner product in feature space has an equivalent kernel in input space, ( , ') ( ), ( ')

K x x x x     (16)
provided certain conditions hold. If K is a symmetric positive definite function, which satisfies Mercer's Conditions, ( , ') ( ), ( '),

m m m m K x x a x x      0 m a  (17) ( , ') ( ) ( ') ' 0, K x x g x g x dxdx   2 , gL  (18)
then the kernel represents a legitimate inner product in feature space. Valid functions that satisfy Mercer's conditions are now given, which unless stated are valid for all real x and x'. A polynomial mapping is a popular method for non-linear modeling as,

Polynomial (homogeneous): ( , ') ( . ') d k x x x x  (19) Polynomial (inhomogeneous): ( , ') ( . ' 1) d k x x x x  (20)
Radial basis functions have received significant attention, most commonly with a Gaussian form:

2 2 ' ( , ') exp( ) 2 xx k x x    (21)
Classical techniques utilizing radial basis functions employ some method of determining a subset of centers. Typically a method of clustering is first employed to select a subset of centers. An attractive feature of the SVM is that this selection is implicit, with each support vector contributing one local Gaussian functions, centered at that data point.

Genetic Algorithm (GA)

This method was inspired by work first reported by Siedlecki and Sklansky [START_REF] Siedlecki | A Note on Genetic Algorithms for Large-Scale Feature Selection[END_REF]. The genetic algorithm is a method for solving both constrained and unconstrained optimization problems that is based on natural selection, the process that drive s biological evolution. In recent years, the genetic algorithms are used in many medical applications [START_REF] Kilic | Colonic Polyp Detection in CT Colonography with Fuzzy Rule based Template Maching[END_REF].

The genetic algorithm repeatedly modifies a population of individual solutions. At each step, the genetic algorithm selects individuals at random from the current population to be parents and uses them to produce the children for the next generation. Over successive generations, the population "evolves" toward an optimal solution. Genetic algorithms can be applied to solve a variety of optimization problems that are not well suited for standard optimization algorithms, including problems in which the objective function is discontinuous, non differentiable, stochastic, or highly nonlinear.

The genetic algorithm uses three main types of rules at each step to create the next generation from the current population. Firstly, selection rules select the individuals, called parents that contribute to the population at the next generation. Secondly, crossover rules combine two parents to form children for the next generation. Finally, mutation rules apply random changes to individual parents to form children. GA is designed to simulate the evolutionary processes that occur in nature [START_REF] Peña-Reyes | Evolutionary Computation in Medicine: an Overview[END_REF][START_REF] Koza | Evolving Inventio ns[END_REF].

The basic idea is derived from the Darwinian theory of survival of the fittest. Three fundamental mechanisms drive the evolutionary process: selection, crossover and mutation within chromosomes. As in nature, each mechanism occurs with a certain probability, allowing for some randomness. Selection occurs on the current population by choosing the most-fit individuals to reproduce. Reproduction, then, can result in the crossover and/or mutation of parent genes to form new solutions. The ratio of heuristic and stochastic decisions, at best, creates a natural balance between survival and evolution.

For each generation of the GA, individual solutions are evaluated using a fitness function.

The evaluation method is a crucial component of the selection process since offspring for the next generation are determined by the fitness values of the present population. Figure 5 provides a simple diagram of the iterative nature of genetic algorithms. The generational process ends when the user-defined goal is reached. In most cases, the number of generations is a constant and is set by the user [START_REF] Miller | Feature Selection for Computer-Aided Polyp Detection using Genetic Algorithm[END_REF]. The aim is to design a decision rule which is easy to compute and yields the smallest possible probability of misclassification of input patterns from the feature space.

Image Processing and Mechanical Assessment

The next section details the techniques developed to investigate bone micro architecture and its mechanical stiffness.

Hybrid Skeleton Graph Analysis (HSGA)

The Hybrid Skeleton Graph Analysis (HSGA) [START_REF] Aufort | Hybrid skeleton graph analysis of disordered porous media. Application to trabecular bone[END_REF] model is a high-level representation of a hybrid shaped porous medium composed of rods and plates such as trabecular bone. The morphology and topology of disordered porous media were first investigated using global methods based on physical models that cannot give complete information about the trabeculae's

Genetic Learning and Evolution

Classifier Feature Extractor structure and its local properties. In 2000 a new method called Line Skeleton Graph Analysis (LSGA) [START_REF] Aufort | A New 3D Shape-Dependent Skeletonization Method[END_REF] was introduced for studying porous media using a skeleton-based technique.

However, as it uses a curve skeleton, the LSGA has its drawbacks since all non-cylindrical shapes are better described by 2D-surfaces rather than by 1D-curves.

a) Hybrid skeleton

The HSGA [START_REF] Aufort | Hybrid skeleton graph analysis of disordered porous media. Application to trabecular bone[END_REF] relies on a new hybrid skeletonization process which is computed by processing curve or surface thinning, depending on the local shape of the object. To switch between the 2 skeleton variants, the HSGA improves a recent algorithm [/***compléter bonnassie] which classifies the voxels of an object according to their topological predisposition to belong to a plate or to a rod zone [START_REF] Aufort | Shape classification techniques for discrete 3D porous media. Application to trabecular bone[END_REF]. Then, surface and curve thinning algorithms are applied depending on whether each primitive belongs to the set of rods or plates. The proposed hybrid thinning algorithm preserves connectivity which is an essential feature when characterizing porous media such as trabecular bone. Figure 8

b) Classification

Once the hybrid skeleton (figure 8(a)) has been computed, a classification step is applied to label each voxel of the skeleton according to its structural role. The classification consists in affecting a class identifier to each voxel of the skeleton according to its structural role. We define 4 classes of voxels: "rod", "line-end", "plate" and "node". As in [START_REF] Aufort | A New 3D Shape-Dependent Skeletonization Method[END_REF], 2 voxels of the solid phase are neighbors if they are 26-connex (i.e. they share at least one corner). Similarly, voxels of the pore phase are neighbors if they are 6-connex (i.e. they share at least a face).

The voxel classification is based on a 3-step process:

1

Step 1: Initialization All voxels are set to "plate".

2

Step 2: Determination of rod-shaped 1D-curves A voxel is marked as "rod" if it has only 2 solid neighbors that are not neighbors themselves. A voxel is marked as "line-end" if it has only 1 neighbor.

3

Step 3: Determination of interfaces between elements Any voxel is set to node if:

-It has more than 2 "rod" or "line-end" neighbors -It is a "plate" and has a "rod" or "line-end" neighbor Figure 8(c) shows the result of the classification step. Each voxel of the skeleton has been marked as belonging to a "rod", "line-end", "plate" or "node" class.

c) Individualization

After classification, the role of each voxel in the skeleton can be determined. However, the plates and rods of the structure cannot be processed one by one, since they have not been extracted and individualized. In each class, the structure elements have to be marked separately. To do so, all the information associated to one plate or rod must be collected. This is the object of the individualization process which processes recursively all the voxels of the structure and determines the interfaces of each element of the object. The elements individualization algorithm is implemented as a 3-step loop:

4

Step 1: Finding a solid phase element An initial voxel (the seed) is chosen in the classified skeleton. It can be of type "plate" for a plate element or "rod" for a rod element.

5

Step 2: Spreading the information

The seed is recursively spread through the skeleton voxels until the element boundaries are found. Boundary voxels can be of type "node" or "line-end".

6

Step 3: Registering the new element

Once an element has been entirely defined, its voxels and interfaces are registered in the model.

this information is registered in the HS GA model. Figure 8(c) illustrates how different elements of the structure are individualized.

d) Segmentation

The classified and individualized skeleton supports a standalone use as it already contains information about the medium. The analysis can be improved however by processing segmentation of the original object. For this purpose, we implement a 3D region-growth iterative process. It takes the skeleton voxels of each element as a seed and iteratively merges neighbors from the original volume. As a result, each solid voxel of the object is finally associated to an element of the HSGA model (plate or rod). The HSGA is used as a basis for the generation of finite element models.

Mechanical Assessment

The reference for Finite Elements (FE) analysis of discrete sa mples is unquestionably voxel-toelement conversion as discussed in many articles [START_REF] Ulrich | Finite element analysis of trabecular bone structure: a comparison of image-based meshing techniques[END_REF][START_REF] Bayraktar | Nonlinear micro finite element analysis of human trabecular bone[END_REF]. However, in the case of skeleton-based models, other types of elements can be used to simplify large-scale problems and prevent time resource consumption. In this section, we explain our modeling choices, relative to the compromise between simplification and loss of accuracy in the simulation. First, the method used to convert rod shapes to beam chains is explained. Then the triangulation technique for converting plate shapes to shell elements is described.

a) Rods to beam elements chains conversion

The FE that matches the geometry of a straight rod is the beam element. Its geometry is described as a 1D segment, which is assigned a circular cross section. This technique was first investigated in [START_REF] Pothuaud | A new computational efficient approach for trabecular bone analysis using beam models generated with skeletonized graph technique[END_REF] to assess the stiffness of trabecular bone. However, results show that modeling the bone by a simple rod network is not sufficient to obtain an accurate assessment of stiffness, due to geometrical lack of accuracy. Lenthe et al. [START_REF] Lenthe | Specimen specific beam models for fast and accurate prediction of human trabecular bone mechanica l properties[END_REF] explored this field of beam-modeling rod shapes in porous media, but did not resolve the shape accuracy issues.

In order to convert a rod item to FE, we introduced the "beam chains" concept [START_REF] Aufort | Nouvelle approche de modélisation de milieux poreux. Application à l'os trabéculaire[END_REF]. Based on a feature extraction technique used in the field of 3D animation [START_REF] Reinders | Skeleton graph generation for feature shape description[END_REF], the beam chain introduces evenly set intermediate nodes on the curve skeleton as seen in figure 9. This process is called "splitting" as it breaks the curve into small segments that better match the curvature of the

b) Plates to shell elements conversion

Plate zones are badly described in the case of beam-only models [START_REF] Aufort | Nouvelle approche de modélisation de milieux poreux. Application à l'os trabéculaire[END_REF][START_REF] Aufort | Shape classification techniques for discrete 3D porous media. Application to trabecular bone[END_REF], which lead to a nonnegligible bias for morphological results [START_REF] Pothuaud | A new computational efficient approach for trabecular bone analysis using beam models generated with skeletonized graph technique[END_REF]. It is suspected that this lack of geometrical accuracy also alters mechanical results. In fact, modeling plates is much more challenging than modeling rods, since research in this area is scanty. Recently, Lenthe et al. [START_REF] Lenthe | Specimen specific beam models for fast and accurate prediction of human trabecular bone mechanica l properties[END_REF] proposed converting a plate into a set of beams instead of a single beam. Ye t, the efficiency of this conversion is questionable since the notion of beam set is not geometrically obvious, and the beam sections are model-dependent. We present here an original approach that combines the power of a new triangulation method and a better choice of FE type to improve plate modeling.

The FE used to describe planar shapes is the shell element. It is defined as a 2D medial surface geometry (either triangle or quad strips for example), on which each element is assigned a thickness value. The challenge in this case was to be able to transform the 26-connected plate voxels sets into a list of simple 2D primitives such as triangles or quads. Meshing iso-surfaces has been widely explored in computer graphics since the early 80s, leading to a great range of techniques: surface fitting, surface tracking (also known as continuation methods) and spatial sampling. However, the case of crossing and stacking manifold surfaces from disordered data is still difficult to handle with criteria like Delaunay tracking, and none of these techniques suited our needs. We chose therefore to draw on a well-known spatial sampling method: Marching Cubes (MC) [START_REF] Lorensen | Marching cubes: a high resolution 3D surface construction algorithm[END_REF], which subdivides space into cells and searches those that intersect the implicit surface. Our algorithm, called Surface Marching Cubes (SMC), computes the full-resolution triangulation of any 26-connected surface by following the MC principle while using a new set of neighborhood patterns. All the triangles linking the voxels of the surface are generated. Figure 10 presents the result of the triangulation of a simple 26-connex surface set. Finally, each shell element (i.e. each triangle) is assigned a section according to the local thickness value in the thickness map, as was done for the beam elements in the c ase of rod shapes. The HSGA is used as a basis for the generation of FE models. Thickness map matching [START_REF] Hilderbrand | A new method for the model independent assessment of thickness in threedimensional images[END_REF] is processed to improve the geometrical accuracy of all the FE models generated, using local cross-section values.

All the FE models were imported into the commercial software "Abaqus", which was used to estimate their apparent Young's Modulus. As the simulation required material definition, we assigned the entire model the same behavior, with a bone tissue characterized by an arbitrary Young's modulus of 15 GPa and a Poisson's coefficient of 0.3. These values are close to those found in the literature [START_REF] Lenthe | Specimen specific beam models for fast and accurate prediction of human trabecular bone mechanica l properties[END_REF][START_REF] Bourne | Finite element models predict cancellous apparent modulus when tissue modulus is scaled from specimen ct-attenuation[END_REF]. In this kind of comparative study, the material's behavior does not really matter since the measured reaction forces are compared relatively. Once defined, each model was tested in compression in the 3 space directions (x, y and z). Boundary conditions on the cube's faces were set to zero for translations perpendicular to the faces and their complementary rotations. Complementary translations were also blocked for the 2 faces perpendicular to the compression axis. Then, a displacement value Δl was applied on the compression direction's front face. A small displacement of 2% of the cube's size (i.e. 0.0913 mm) avoided non-linearity issues. The apparent Young's modulus (Eapp) of the model was then computed [START_REF] Lenthe | Specimen specific beam models for fast and accurate prediction of human trabecular bone mechanica l properties[END_REF] as:

2 ( / ) /( / ) Eapp RF l l l    (22)
where ΣRF represents the measured sum of the Reaction Forces (RF) on each node the compression face, and l is the size of the cube's side.

Results

The HSGA is an efficient tool for the analysis of porous media. This section presents its usefulness to describe trabecular bone samples from a clinical study. For this purpose, the medical staff at the hospital of Orleans (France) provided us with 2 populations of 9 sa mples each, extracted from post-mortem femoral head and acquired using a Sky Scan micro-scanner with a high-resolution µ-CT (12 µm). The first 9 samples came from OsteoArthritic (OA) patients, whose bone structure is known to be hypertrophied, increasing the bone density. The other 9 samples were extracted from OsteoPorotic (OP) patients characterized by the deterioration of bone micro architecture which leads to bone fragility and fracture risk. As the characteristics of the 2 populations are previously known, they can be used to verify the separating power of any feature that is said to reflect the bone micro architecture alterations. The numerical samples are 400 3 isotropic voxels 8-bit grey level volumes. They were pre-processed as follows: after applying a median filter, they were binarized at the local minimum threshold between the 2 modes of their histogram. The Hoshen-Kopelman [START_REF] Hoshen | Percolation and cluster distribution i. cluster multiple labeling technique and critical concentration algorithm[END_REF] clustering algorithm was then used to remove non-connected solid voxel sets, as there can be no isolated material in the bone sample. Figure 11 shows 2 extracts from an OA and an OP sample. It can be seen that there is more solid material in the OA sample.

(OA) (OP) Figure 12 shows that none of the computed 27 parameters is able to distinguish significantly between the 2 populations studied. This study involves 18 samples with 27 input parameters and 1 output parameter. 12 samples are used for prediction and 6 for the te st to be classified as OA or OP with Adaptive Neuro-Fuzzy Inference System (ANFIS), Support Vector Machines (SVM)

and Genetic Algorithm (GA). Firstly we used ANFIS (Figure 13).

Figure13. Architecture of our ANFIS model: 27 inputs and 1 output with 4 rules.

In ANFIS model, Gaussian Membership Function and 5 rules were used. Therefore, 18

samples are not enough for accurate prediction with ANFIS. We could not obtain satisfactory results and it is shown in figure 14. 18 samples are not enough for accurate prediction with ANFIS. We could not obtain satisfactory results using ANFIS, as shown in figure 14. We applied the same real data group for SVM with Linear kernel function. For linear kernel C=1, Epsilon = e -12 and tolerance parameter = 0.100 are chosen. For Genetic Algorithm; population size is 50, crossover rate is 0.9 and mutation rate is 0.07.

Classification of OA and OP using Genetic Algorithm and SVM Algorithm are simulated for various training-test ratios (%) of all data set such as 50-50 and 66.6-33.3. In order to obtain better network generalization, 4-fold cross validations are also used.

Discussion

This paper presents an original association of image processing and artificial intelligence methods for a better classification of two studied populations composed of 9 osteoarthritis and 9

osteoporosis samples denoted OA and OP respectively.

Concerning image processing, we have presented a precise assessment technique for the 3D characterization of a complex porous medium such as trabecular bone. Hybrid Skeleton Graph Analysis (HSGA) is a technique which generates structural models of the object. It consists in an efficient combination of curve and surface thinning techniques and considers the local shape of each element that composes the structure of the medium. HSGA can classify each element of the structure as rod-like or plate-like without ignoring the intersections and the termini. Features extracted from the HSGA contain significant topological and morphological information. The HSGA is used to extract mechanical features thanks to a finite element analysis.

Twenty-seven features are processed on each of the 18 bone samples. A comparative discriminate analysis using artificial intelligence is applied to identify each of the samples as OA or OP.

For the first test, Adaptive Neuro-Fuzzy Inference System (ANFIS) was used but its success in identifying each of the samples was poor due to lack of samples and many nodes (Figure 14). We then applied Support Vector Machines (SVM) with Linear Polynomial Kernel and obtained a high degree of success. Finally, we used Genetic Algorithm (GA) and obtained correct predictions and 100 % success (Table 2).

The results presented in this paper are consistent with our expectations, since geometrical improvement in the models leads to better and more precise characterization of the structure. The long-term aim of our work is to develop a biomechanical simulation protocol that could be used to virtually characterize the stiffness of trabecular bone samples [START_REF] Aufort | Mechanical assessment of porous media using hybrid skeleton graph analysis and finite elements. Application to trabecular bone[END_REF]. This will contribute significantly to detecting bone fragilities simply by acquiring 3D images, for example with a high-resolution in-vivo CT and simulating mechanical compression using fast and precise FE models. Associating different techniques with artificial intelligence could be used by physicians to complete their diagnosis efficiently. The drawback of this study is the small number o f samples in the populations. Further work is to be carried out on larger sets of trabecular bone samples.

Figure 1 .

 1 Figure 1. ANFIS architecture.

c

  represents the membership function's (Mfs) center and  determined the Mfs width.

Figure 2 .

 2 Figure 2. Mapping the Input and Output space.

Figure 3 .

 3 Figure 3. General separation of optimal hyper plane.

Figure 4 .

 4 Figure 4. Linear classifier and margins.

  is a kernel function which high-dimensional feature space.

Figure 5 .

 5 Figure 5. Flow diagram depicting the evolutionary process of genetic algorithm.

Figure 6 .

 6 Figure 6. General separation of two classes by GA.

Figure 7 .

 7 Figure 7. Feature extractor and classifier with feedback learning system.

  (a) and figure8(b) show respectively an extracted portion of a trabecular bone sample and its hybrid skeleton.

Figure 8 (Figure 8 .

 88 Figure 8. Illustration of the HSGA model and its different computing steps on a trabecular bone sample. Original trabecular bone (a), hybrid skeleton (b), classified hybrid skeleton (c), final HSGA segmented model (d). In (b) and (c), the skeleton is stacked over the original object.

Figure 9 .

 9 Figure 9. Illustration of the rods to beam elements chains conversion: voxels of the 1D path extracted from its curve skeleton with intermediate nodes (a), simple cylinder model assumption using the rod's volume and skeleton length to compute its section as S = V / l (b), and beam chain modeling of the item using local thicknesses (c)

Figure 10 .

 10 Figure 10. Illustration of the plates to shell elements conversion: 26-connex voxels of the 2D surface extracted from its surface skeleton (a), and example of a shell elements triangulation using Surface Marching Cubes (b).

Figure 11 .

 11 Figure 11. Two 64 3 voxels extracts from Osteoarthritis (OA) and Osteoporotic (OP) samples to illustrate the micro architectural differences in the two trabecular bones. The solid phase (bone) is illustrated in black while the pore phase is in white.

Figure 12 .

 12 Figure 12. Input and output parameters. a) Distribution of 27 Input parameters versus 1 Output Parameter, b) One output Parameter of two classes (OA, OP).

Figure 14 .

 14 Figure 14. Prediction after training, "o" is the actual value and "*" is the predicted value.

Table 1 :

 1 Input parameters obtained from HSGA and FE Analysis.

	b) Morphological parameters
	From a morphological point of view, the HSGA enables the following non exhaustive list of
	parameters to be measured: Bone Mineral Density (BMD, defined as Bone Volume over Total
	Volume, BV/TV), Bone Surface over Total Volume (BS/TV, mm -1 ), Bone Surface over Bone
	Volume (BS/BV, mm -1 ), Rod Volume (Rod.V, mm 3 ), Rod Proportion (Rod.Prop, % ), Plate
	Volume (Plate.V, mm 3 ), Plate Proportion (Plate.Prop, %), Element Number (El.N defined as
	Ro.N plus Pl.N), Plate Number (Pl.N), Rod Number (Ro.N), Rod Length (Ro.L, mm), Rod
	Volume (Ro.V, mm 3 ), Rod Section (Ro.S, mm 2 ), Rod Thickness (Ro.Th, mm), Plate Surface
	(Pl.S, mm 2 ), Plate Volume (Pl.V, mm 3), Plate Thickness (Pl.Th, mm).

c) Mechanical parameters

Finally, associating the HSGA model to Finite Element Analysis allows measurement of new features such as: Mesh Shell Number (M.Sh.N), Mesh Node Number (M.No.N), Mesh Beam Number (M.Be.N), and Young's Modulus (E App (x, y, z)) in the 3 Space directions (x,x and z).

Table

1

gives all the input parameters with minimum, maximum, mean and standard deviation values for each of the 27 parameters described above. Figure

12

shows all the input parameters' distribution versus one output parameter value.

Table 2 .

 2 Table 2 lists confusion matrix and obtained classification success rate. Confusion Matrix and results of simulated methods.

	Methods		Evaluate on 9 training and 9 testing data		4-fold cross-validation	Evaluate on 12 training and 6 testing data
		OA OP Success Rate=100%	OA	OP	Success Rate=88.88%	OA	OP	Success Rate=83.33%
	SVM	4	0	OA	8	1	OA	3	0	OA
		0	5	OP	1	8	OP	1	2	OP
		OA OP Success Rate=100%	OA	OP	Success Rate=94.44%	OA	OP	Success Rate=100%
	GA	4	0	OA	9	0	OA	3	0	OA
		0	5	OP	1	8	OP	0	3	OP


Rachid Jennane
email: rachid.jennane@univ-orleans.fr

Rachid Harba

Gérald Lemineur

Stéphanie Bretteil

Anne Estrade

Claude Laurent Benhamou

Estimation of the 3D self-similarity parameter of trabecular bone from its 2D projection

Keywords: Self-similarity, fractional Brownian motion, fractal, trabecular bone, osteoporosis, X-ray images, computed tomography

It has been shown that the analysis of two dimensional (2D) bone X-ray images based on the fractional Brownian motion (fBm) model is a good indicator for quantifying alterations in the three dimensional (3D) bone micro-architecture. However, this 2D measurement is not a direct assessment of the 3D bone properties. In this paper, we first show that S 3D , the self-similarity parameter of 3D fBm, is linked to S 2D , that of its 2D projection, by S 3D = S 2D -0.5. In the light of this theoretical result, we have experimentally examined whether this relation holds for trabecular bone. Twenty one specimens of trabecular bone were derived from frozen human femoral heads. They were digitized using a high resolution -CT. Their projections were simulated numerically by summing the data in the three orthogonal directions and both 3D and 2D self-similarity parameters were measured. Results show that the self-similarity of the 3D bone volumes and that of their projections are linked by the previous equation. This demonstrates that a simple projection provides 3D information about the bone structure. This information can be a valuable adjunct to the bone mineral density for the early diagnosis of osteoporosis.

Introduction

Osteoporosis is characterized by low bone mass and micro-architectural deteriorations of bone tissue [START_REF] Parfitt | Implications of Architecture for the Pathogenesis and Prevention of Vertebral fractural[END_REF]. These modifications lead to enhanced bone fragility, and consequently increase the risk of fracture (Consensus Development Conference, 1991). The distribution of mineral content in three dimensions, and the loading of the skeleton are additional factors that contribute to bone fragility.

The diagnosis of osteoporosis is mainly based on dual energy X-ray absorptiometry which amounts to measuring bone mass. Characterization of the trabecular structural properties appears to be an important adjunct to the measurement of bone mass in determining fracture risk with greater accuracy [START_REF] Goldstein | The Mechanical Properties of Trabecular Bone: Dependence on Anatomic location and Function[END_REF]. Using histological and stereological analysis, it has been shown that, by combining structural features with bone density, nearly all of the variability in mechanically measured Young's moduli could be explained [START_REF] Goulet | The relation between the structural and orthogonal compressive properties of trabecular bone[END_REF][START_REF] Croucher | Assessment of cancellous bone structure: comparison of strut analysis, trabecular bone pattern factor, and narrow space star volume[END_REF]. However, the evaluation of bone structure, by non-invasive procedures, remains a difficult issue [START_REF] Faulkner | Noninvasive Measurements of Bone Mass, Structure and Strength: Current Methods and Experimental Techniques[END_REF]. Bone X-ray radiographs have been suggested as a means of quantifying trabecular changes (Rockoff et al., 1986;[START_REF] Aggarwal | A survey of osteoporosis using the calcaneum as an index[END_REF]. This method shows potential as we know that radiographs are feasible and are not life threatening. Among many methods that can extract the information from natural images, fractals are efficient candidates [START_REF] Pentland | Fractal-Based Description of Natural Scenes[END_REF], and fractional Brownian motions (fBm) [START_REF] Mandelbrot | Fractional Brownian Motion, Fractional Noise and Applications[END_REF]) are widely used. FBm is governed by the single Hurst parameter, H (0 < H < 1). It has been shown that H and the intuitive notion of roughness are closely linked: the lower H is, the higher is the roughness [START_REF] Pentland | Fractal-Based Description of Natural Scenes[END_REF]. Concerning bone radiographs, this idea was used in [START_REF] Lundhal | Fractional Brownian motion : a maximum likelihood estimator and its application to image Texture[END_REF] to characterize the texture roughness. It was shown that H follows the clinically expected course of osteoporosis closely, in that a lower H indicates a greater degree of osteoporosis. Indeed, during osteoporosis, trabeculae become thinner, perforations occur, and some of the fine trabecular structures are broken. Radiographs appear naturally rougher, leading to a decrease in H.

These results have been confirmed by subsequent studies [START_REF] Benhamou | Fractal Organization of the Trabecular Bone Images on Calcaneus Radiographs[END_REF][START_REF] Majumdar | Fractal Analysis of Radiographs: Comparison of Techniques and Correlation with BMD and Biomechanics[END_REF][START_REF] Majumdar | In vivo assessment of trabecular bone structure using fractal analysis of distal radius radiographs[END_REF][START_REF] Prouteau | Fractal analysis of bone texture: a screening tool for stress fracture risk?[END_REF]. Clinical results have also pointed out the strong potential of the H parameter of fBm in quantifying osteoporosis on bone X-ray radiographs [START_REF] Pothuaud | Fractal Analysis of trabecular Bone Texture on Radiographs: Discriminant value in post menopausal Osteoporosis[END_REF][START_REF] Lespessailles | Biomechanical Properties of human os calcanei : relationships with bone density and fractal evaluation of bone microarchitecture[END_REF].

An X-ray radiograph is a two-dimensional (2D) projection of a three-dimensional (3D) structure.

The difficulty in using such 2D projectio ns is in determining exactly how the 2D parameters measured on radiographs are related to 3D ones. This problem is part of the more general one of using 2D images, either slices or projections, to characterize 3D materials. This is a difficult issue, though some results are available in the case of fractals. For example, it is well known that a 2D slice of a fractal volume has the same H parameter as the volume itself [START_REF] Mandelbrot | The fractal geometry of nature[END_REF]. However, for projections, no equivalent link has been established. One way to elucidate the problem is to conduct experimental studies on real material [START_REF] Pothuaud | Fractal Dimension of Trabecular Bone Projection Texture Is Related to Three-Dimensional Microarchitecture[END_REF]Jennane et al., 2001;[START_REF] Guggenbuhl | Texture analysis of X-ray radiographs of iliac bone is correlated with bone micro-CT[END_REF]. However, this approach which consists in estimating 3D properties from 2D projections suffers from the lack of a strong theoretical background. Moreover, the relations between 2D and 3D parameters remain unclear.

Recently, an isotropic analysis for Gaussian models having stationary increments has been suggested [START_REF] Bonami | Anisotropic Analysis of some Gaussian Models[END_REF]. This study was applied in (Jennane et al., 2001) where we proposed an oriented fractal analysis for anisotropic fractional Brownian motion. The theoretical study presented in [START_REF] Bonami | Anisotropic Analysis of some Gaussian Models[END_REF] can also be applied to 3D bone images. In this paper, we first establish a formal link between S 3D , the self-similarity of a 3D continuous fBm process and S 2D , the self-similarity of its 2D projection, that is S 2D = S 3D + 0.5. In the light of this theoretical result, we experimentally examine whether or not this formal link holds for trabecular bone. With this aim in view, we study high resolution micro computed tomography (-CT) images and simulate their 2D projections obtained by summation of the voxels of the 3D structure. Then, 2D and 3D self-similarity parameters are measured and compared. If the formal link S 2D = S 3D + 0.5 holds for trabecular bone, information on 3D structure can be extracted out of 2D radiographs and can help in the early diagnosis of osteoporosis.

The paper is organized as follows: In section II, we de monstrate the relation between 3D and 2D self-similarity parameters for fBm. Section III presents the preparation and acquisition of the μ -CT bone images. The measurement methods for 3D and 2D self-similarity descriptors are explained in section IV. Section V presents the results and we conclude with a brief discussion of this work.

Theory

It has been shown that the H parameter of the fBm is a good indicator for the alterations of the 3D bone micro-architecture seen on 2D bone radiographs [START_REF] Benhamou | Fractal Organization of the Trabecular Bone Images on Calcaneus Radiographs[END_REF][START_REF] Majumdar | Fractal Analysis of Radiographs: Comparison of Techniques and Correlation with BMD and Biomechanics[END_REF].

However, this 2D measurement is not a direct assessment of the 3D bone properties. Based on a recent paper [START_REF] Bonami | Anisotropic Analysis of some Gaussian Models[END_REF], it is possible to establish a formal link between S 3D , the selfsimilarity of a 3D continuous fBm and S 2D , the self-similarity of its 2D projection. The aim of this section is to demonstrate this relation.

Let us first recall a result and the way to obtain it. We are interested in a signal indexed by the position in the 3D space domain, given by a fBm, B H : (B H (x); x in R 3 ). One of the characteristics of fBm is self-similarity, that is to say:

for all scales  (non-negative),     3 3 ; ) ( ; ) (      x x B x x B H H H   , (1)
where  means equal in distribution. The Hurst parameter H is also the 3D self-similarity index denoted, S 3D . Thus, S 3D = H.

The self-similarity identity (1) can be obtained as follows. Since B H is a zero mean Gaussian field with stationary increments, its distribution is completely defined by the variance of its increments

B H (x): 3 2 in x all for , )) ( (    H H x x B Var , (2)
where, for a multidimensional lag x at any point x 0 in  3 ,

). ( ) ( ) ( 0 0 x B x x B x B H H H    
(3)

x denotes the Euclidian norm in  3 : 2 3 2 2 2 1 3 2 1 for x x x x ) , x , , x (x x     . (4)
Now, introducing a change of scale xx gives the following:

)

) ( ( )) ( ( 2 2 2 2 x B Var x x x B Var H H H H H H          , ( 5
)
which shows (1).

In the following, the term projection will be used for a simple sum of the planes of the 3D structure.

To establish the formal link between the self-similarity of a 3D fBm and that of its 2D projection, we will now consider the projection along one of the three orthogonal directions. This direction plays a special role and we will denote for clarity x = (x 1 , x 2 , z), x   3 the 3D position in the space domain where z is the projection direction. Identically,  = ( 1 , 2,  z ),    3 is the 3D position in the frequency domain. When B H is projected along the z axis, one obtains a 2D Gaussian field Y H :

   dz z z x x B x x Y H H ) ( ) , , ( ) , ( 2 1 2 1  . (6)
(z) is a window that ensures the convergence of the previous integral. Y H is no longer a fBm but still has stationary increments as proved in the following. Consider the spectral representation of a 3D fBm [START_REF] Reed | Spectral Representation of Fractional Brownian Motion in n Dimensions and its Properties[END_REF]:

           3 2 2 1 1 , , ) , , ( 1 1 ) , , ( 2 1 2 3 2 1 ) ( 2 1 z H z z x x i H H dW e C z x x B z          , ( 7
)
where C H is a non negative constant and dW is a 3D Brownian measure in frequency. The projected field can be represented as follows:

             3 2 2 1 1 , , ) ( ) , , ( 1 1 ) 0 , 0 ( ) , ( ) , ( 2 1 2 3 2 1 ) ( 2 1 2 1 z z H z x x i H H H H dW e C Y x x Y x x Y           , ( 8
)
where ˆ is the Fourier transform of . Computing the variance of Y H we obtain

  . 1 ) ( ˆ ) ( 2 1 sin 1 )) , ( ( 2 2 1 2 3 2 2 2 2 1 2 2 2 1 1 2 2 H 2 1                                     d d d x x C x x Y Var z H z z H (9)
This variance does not satisfy any identity such as (5) and therefore the projected field Y H is not selfsimilar. However, we are going to demonstrate that Y H is self-similar for small scales.

Let us determine the variance of Y H (x 1 , x 2 ). Introducing the change of variable  

  . 1 ) ( ˆ ) ( 2 1 sin )) , ( ( 2 2 1 2 3 2 2 2 2 2 1 2 2 2 1 1 2 2 1 2H 2 1                                          d d d x x C x x Y Var z H z z H H (10)
For small scale , i.e.   0, and if H < 0.5:

    H 0.5 H 2 2 1 2 2 3 2 2 2 1 2 2 1 1 2 2 H 1 2H 2 1 ) ( ˆ 1 ) ( 2 1 sin C )) , ( ( 2 V d d d x x x x Y Var z z H H                                        (11)
where V H is a constant which does not depend on . In that case, by analogy with identity (5), we say that Y H is asymptotically self-similar of index H + 0.5. In other words, S 2D = H + 0.5. Actually, the previous result is only valid for H < 0.5. For H > 0.5, the result still remains true when using the 2-order increments of Y H , that is to say

H 1 2H 2 V )) ( (      x Y Var H , (12)
where

), ( ) ( 2 ) 2 ( ) ( 2 x Y x x Y x x Y x Y o o       (13)
with V' H a constant which does not depend on . The latter equation is valid for any point x o in  2 .

Hence, the asymptotical self-similarity index of the projection Y H which is a 2D signal process is S 2D = H + 0.5 while the self-similarity index of the original 3D process is S 3D = H. For small scale , we can write that:

0.5. - S S 2D 3D  (14)
Several remarks can be made at this point.

Remark 1: Although this result has been shown in the 3D case, it is valid for any space dimension above one.

Remark 2: The increase in the self-similarity of the projection is not unexpected since this operation is an integration process. This action leads to softer objects, thus increasing the self-similarity index.

Remark 3: For the main result ( 14) to be valid, the self-similarity parameter of the projection should be estimated for small scales , or in other words for high frequency. The practical estimation of S 2D should take this fact into consideration.

Remark 4: The projected process is not self-similar but asymptotically self-similar at high frequency of index H + 0.5 in the interval ]0.5, 1.5[. For the large scale domain, i.e. the low frequency domain, it ca n be shown that the projected process Y H is asymptotically self-similar of parameter H [START_REF] Bonami | Anisotropic Analysis of some Gaussian Models[END_REF]. As a conclusion, Y H has 2 different asymptotical self-similar parameters, H at low frequency and H + 0.5 at high frequency.

We will now explore if this theoretical result holds for 3D trabecular bones and their 2D projection.

The next paragraph explains how the bone samples were acquired and digitized.

Preparation and acquisition of bone specimens

A total of 21 frozen human femoral heads were used to derive 21 specimens of trabecular bone.

Cylindrical samples were prepared under continuous water irrigation using a precision diamond circular saw. The samples were oriented with respect to anatomic axes and cut to a dimension of 6 mm thickness and 8 mm diameter. All the samples were defatted chemically (several cycles of submerging in dichloromethane).

Images were obtained using the Skyscan 1072 high-resolution -CT. The X-ray source was set at 80kV and 100µA, and the magnification was fixed to get a pixel size of 12 µm. A 10241024 12-bit digital cooled CCD coupled to a scintillator was used to record the radiographic projections. 209 projections were acquired over an angular range of 180° (angular step of 0.9°). Due to the cone beam the radiographic images were processed with the Feldkamp algorithm. The exposure time for one radiographic projection was about 5.9s so that the total scan for each sample lasted 1 hour. All the radiographic images were transferred to a personal computer (PIII, 933MHz) and the image slices were reconstructed using the conebeam reconstruction software version 2.6.

We used the central part of each image which is 4.80 mm 3 , equivalent to 400 3 pixels. Fig. 1 shows two such cubes having different 3D structures. The differences between these two samples are obvious: low density and connectivity for Fig. 1a, high density and connectivity for Fig. 1b. This is confirmed by the following:

-BV/TV (the ratio between the volume of trabecular bone (BV) and the total volume (TV)) was 0.182 for sample (a) and 0.306 for (b) -Conn.D (the connectivity density parameter) was -0.0044 for sample (a) and -0.0075 for (b).

-For the 21 samples, the mean and standard deviation for BV/TV was 0.210  0.066, and -0.0111  0.0096 for Conn.D.

We now wish to compute the projection of these bone volumes. To do so, we used a simple projection which consists of summing the data in the three orthogonal directions x 1 , x 2 and z. Fig. 2 illustrates two projections along the z axis of the two bones presented in Fig. 1. As can be seen, the projection of the bone volumes generates grey level images. It is also noticed that image 2.a is different from image 2.b as is the case for the bone volumes of Fig. 1. It is then expected that the 3D self-similar parameter and that of its 2D projection are linked.

Our aim is to apply to a binary object (trabecular bone) a result concerning the relation between 3D

and 2D self-similarity demonstrated using a continuous fractal model (fBm). A better approach would be to use a 3D binary fractal model indexed by R 3 taking values in {0, 1}. At this point of our project, such work is under active investigation. We hope that in a near future we will be able to propose such a result for some class of binary fractals.

We will now present the chosen tools to access the self-similarity of these 3D and 2D images.

Measure ment of 2D and 3D self-similarity

A wide variety of computer algorithms for estimating the fractal index of a structure exist. The most widely used algorithms in signal and image processing are described in [START_REF] Barnsley | The Science of Fractal Images[END_REF]Jenna ne et al., 1996Jenna ne et al., , 2001)).

Let us now describe the two methods used to compute the 2D and 3D self-similarity parameters.

2D self-similarity: the variance method

The projections of 3D bones are 2D grey level images as seen on Fig. 2. Thus equation ( 5) can be applied to assess the 2D self-similarity. This method is known as the variance method of Pentland [START_REF] Pentland | Fractal-Based Description of Natural Scenes[END_REF]. It requires that the increments I(x) of the grey levels I(x) be computed for different lags . If self-similar of index S 2D , the following relation holds:

)

) ( ( )) ( ( 2 2 x I Var x I Var D S      , ( 15
)
where x is a position index.

In this case, the Log-Log plot of Var(I(x)) versus  is a straight line of slope 2S 2D . Such representations for the variance of the grey level increments obtained over images lines presented in Fig. 2 are illustrated on Fig. 3. Note that  corresponds in the discrete case to the number of pixels considered for the lag's increment. For example, for the 8 first points the spatial resolution varies between 12 μm and 96 μm. It should be remembered that, even if the 3D structure is self-similar, the projected field is only asymptotically self-similar at small lags . This is typically what we observe. For that reason, we choose to estimate S 2D over the first eight points of the Log-Log plot which show a fair linearity. This area characterizes the bone structure at the trabecular thickness scale [START_REF] Harba | Determination of Fractal Scales on Trabecular Bone X-Ray Images[END_REF]. Knowing that the pixel size is about 12 m defines physically a lag whose size is 96 m. This ensures a lag of measure lower than the trabecular size (about 100 to 200 m for the bigger trabeculae) and allows an examination of what is going on inside the trabecular structure to detect perforation or rupture of the trabeculae. S 2D was computed as the mean value of the 400 values estimated over 400 lines of each image.

3D self-similarity: the box counting method

The self-similarity of a continuous object, as for example our 2D images, is easily assessed using equation ( 5). However, concerning our application, it is obvious that 3D trabecular bone is not a continuous object, but rather a binary one. Nevertheless, the notion of self-similarity does exist for such data. Binary fractals are unusual geometric structures that can be used to analyze many biologic structures not amenable to conventional analysis. The full-sized object is composed of elements, these elements composed of the same elements even smaller, and so on down to the microscopic level.

Trabecular bone has a branching pattern, as seen on Fig. 1. One can also see that it exhibits selfsimilarity, that is the trabeculae and the narrow spaces between them look very similar no matter what their size.

The notion of self-similarity is linked to the volume of the structure defined by the following: the original 3D structure is regularly cut in small cubes, or in boxes, of size ρ. The number N(ρ) of boxes containing at least one mineral element is evaluated. The corresponding volume V(ρ) is given by: 3

) ( ) (    N V  (16)
If a binary structure is self-similar of index S 3D then the following relation holds for any scale index

 > 0: ) ( ) ( 3    V V D S  . ( 17
)
Rewriting V(ρ) from ( 16) and replacing V(ρ) in equation ( 17), as a consequence, it easily results that N(ρ) is proportional to

) 3 ( 3D S   
. Returning to fractal objects, the relation between N(ρ) and ρ is (Falconer, 1987):

D C N    ) ( . ( 18
)
where D is the fractal or box dimension of the object and C a constant. Thus, by combining this result with previous ones, the self-similar parameter S 3D is linked to D by:

D S D   3 3 . ( 19
)
This analysis is called the box counting method. The idea that the fractal dimension of trabecular bone might be related to bone architecture and bone strength is an appealing one, since the fractal dimension can be assessed from clinical images. To avoid the phenomenon of length described in [START_REF] Tricot | Evaluation de la dimension fractale d'un graphe[END_REF][START_REF] Feder | Fractals[END_REF], the size of the boxes varies by a power of 2. To be consistent with the measure of S 2D , the slope was measured on the fourth first points of each Log-Log plot which corresponds to a spatial resolution between 12 μm and 96 μm. For both 2D and 3D analysis, the image files were processed with custom designed software developed at our institution using Microsoft Visual C++ language.

Results

Both 2D and 3D self-similarity indices were estimated on the 21 bone volumes and on their projections using the above methods. We examined experimentally if these parameters are linked by S 3D = S 2D -0.5 as is the case for continuous fBm. Fig. 5 shows as an example the evolution of S 2D along the z axis versus S 3D for the 21 bone samples as well as the straight line S 2D = S 3D + 0.5. We notice that S 2D and S 3D are not far from being related by S 2D = S 3D + 0.5. To confirm this tendency, the mean offset, S 2D -S 3D , and its standard deviation are presented for the three orthogonal directions x 1 , x 2 and z in Table 1.

Table 1

Mean and standard deviation for S 2D -S 3D along the three orthogonal directions computed on the 21 bone cubes. This confirms that the above relation is true for trabecular bone. It provides an interesting tool to quantify the 3D architecture from its two dimensional projection. The main advantage is that 2D bone radiographs are not life threatening, cheap and easily obtainable.

Conclusion

In the light of a new theoretical result shown in this paper, this study on high resolution μ-CT bone images has aimed to establish the possible links between the self-similarity of 3D bone volumes and that of their 2D projections. Results show that it is feasible to access the self-similarity of 3D bone volumes solely from their projections computed numerically by summing the original data in the three orthogonal directions. The object of such investigations is a crucial issue to assess changes in bone structure only from 2D X-ray bone radiographs which is cheap and not life threatening. This evaluation of bone micro-architecture may be a useful indicator in a variety of tasks of clinical interest, especially in the prediction of diseases such as osteoporosis. Beyond its application to bone, this result may be interesting for other porous materials. This study opens the way to further theoretical and practical developments.

First, the relation between the three and two-dimensional self-similarity has not been shown for binary fractals. Moreover, it will be interesting to establish a formal link between the self-similarity and the connectivity of the bone which is one of the most useful parameters to assess bone strength. To address these two theoretical issues, we are considering binary self-similar models that can be used for such investigations.

Second, we have shown that the 2D projection of 3D fBm of H parameter in ]0, 1[ is asymptotically self-similar of parameter H at low frequency and asymptotically self-similar of parameter H + 0.5 at high frequency. Hence, the 2D projection of a 3D fBm is a natural example of the generalized model proposed in [START_REF] Kaplan | Extending Self-Similarity for Fractional Brownian Motion[END_REF][START_REF] Perrin | Piecewise fractional Brownian motion[END_REF] taking into account this piecewise behavior.

This theoretical background justifies a local measurement of the H parameter as well as a value of the self-similarity index that takes values above 1. This new approach could be used for a better assessment of the H parameter of bone radiographs.

Finally, from a practical point of view, we are actively developing a high resolution digital X-ray bone prototype at a pixel size of 50 μm which will be used for the clinical evaluation of bone microarchitecture. By associating these theoretical and practical advances, and combining both microarchitecture and density parameters of the bone, we believe that the early and efficient diagnosis of osteoporosis will be achieved within the next few years.

  above, we get:

Fig. 1 .

 1 Fig. 1. two 3D specimens.

Fig. 2 .

 2 Fig. 2. Two projections along the z-axis for the two specimens presented in Fig. 1.

Fig. 3 .

 3 Fig. 3. Log-Log plot of the variance of the increments versus the lag  for the two projections of Fig. 2.

Fig. 4

 4 Fig.4shows a bone volume specimen with d ifferent boxes' edges and the corresponding Log-Log

Fig. 4 .

 4 Fig. 4. (a) A 3D bone volume with several grids placed over the bone structure, and (b) the Log-Log

Fig. 5 .

 5 Fig. 5. Representation of S 2D along the z-axis versus S 3D for the 21 bone samples, as well as the straight

S


Shoko Araki

Francesco Nesta

Emmanuel Vincent

Zbynek Koldovsky

Guido Nolte

Andreas Ziehe

Alexis Benichoux

The 2011 Signal Separation Evaluation Campaign (SiSEC2011): -Audio source separation

. 7 The same algorithm as S2 with additional Binary Masking post-processing. 8 The same algorithm as S4 with additional TF post-processing. 9 The same algorithm as S6 with additional Wiener-Filter like post-processing. 10 The same algorithm as S5 with different parameter settings.

Introduction

The Signal Separation Evaluation Campaign (SiSEC) is a regular campaign focused on the evaluation of methods for signal separation. It was built on the experience of previous evaluation campaigns (e.g., the MLSP'05 Data Analysis Competition1 , the PASCAL Speech Separation Challenge [START_REF] Cooke | Monaural speech separation and recognition challenge[END_REF], and the Stereo Audio Source Separation Evaluation Campaign (SASSEC)) and has been organized since 2008 [START_REF] Vincent | The Signal Separation Evaluation Campaign[END_REF]. SiSEC is not a competition but a community-based scientific evaluation whose aspects are publicly defined. A call for participation precedes the evaluation and aims to define datasets, tasks and evaluation criteria.

This article describes the audio part of SiSEC 2011. In response to the feedback received at SiSEC2008 and SiSEC2010, previous datasets were reorganized as follows:

1. datasets sharing similar scenarios were merged in order to remove some redundancies (e.g. the 2-channel 1-source dataset of the "Source separation in the presence of real-world background noise" task of SiSEC2010 was merged with a new dataset from the PASCAL CHiME Challenge [START_REF] Christensen | The CHiME corpus: a resource and a challenge for computational hearing in multisource environments[END_REF]). 2. tasks with little participation in the previous campaign were excluded; 3. unrealistic data was removed (e.g. the synthetic mixtures of the "Underdetermined speech and music separation" task of SiSEC2010 were eliminated and fresh real-world data was provided).

In general the new campaign was designed so as to better match with real-world scenarios. We believe this data could be of high potential interest for many audio applications in the future years. Specifically, the new datasets embody more realistic features such as a) more reverberant rooms b) real-world diffuse or rapidly varying noise c) source movements. Datasets and tasks are specified in Section 2 and the obtained outcomes are summarized in Section 3. Due to the variety of the submissions, we focus on the general outcomes of the campaign and ask readers to refer to http://sisec.wiki. irisa.fr/ for further details.

Specifications

This section describes the tasks, datasets and evaluation criteria, which were specified in a collaborative fashion. A few initial specifications were first suggested by the organizers. Potential participants were then invited to provide feedback and contribute additional specifications through the wiki or the mailing list.

Tasks

For each dataset, audio mixtures spanning a variety of mixing conditions are provided. The channels x i (t) (1 ≤ i ≤ I) of each mixture signal were generally obtained as x i (t) = ∑ J j=1 s img ij (t), where s img ij (t) is the spatial image of source j (1 ≤ j ≤ J) on channel i [START_REF] Vincent | The Signal Separation Evaluation Campaign[END_REF]. For point sources, s img ij (t) = ∑ τ a ij (t-τ, τ )s j (t-τ ) where s j (t) are the source signals and a ij (t, τ ) the (possibly time-varying) mixing filters. For these mixtures, we specified the following four tasks: T1 Source counting T2 Source signal estimation T3 Source spatial image estimation T4 Source DOA estimation These tasks consist in finding, respectively: (T1) the number of sources J, (T2) the source signals s j (t), (T3) the spatial images s img ij (t) of the sources for all channels i, and (T4) the direction of arrival (DOA) of each source. Participants were asked to submit the results of their systems for T2 and/or T3, and optionally for T1 and/or T4.

Two oracle systems were also considered for benchmarking task T3: ideal binary masking over a short-time Fourier transform (STFT) [START_REF] Vincent | Oracle estimators for the benchmarking of source separation algorithms[END_REF] (O1) and over a cochleagram [START_REF] Wang | On ideal binary mask as the computational goal of auditory scene analysis[END_REF] (O2). These systems require the true source spatial images and provide upper bounds on the performance of binary masking-based systems.

Datasets

Four distinct datasets were provided for SiSEC2011:

D1 Under-determined speech and music mixtures

This dataset includes the stereo dataset D1 from SiSEC2010 [START_REF] Araki | The 2010 signal separation evaluation campaign (SiSEC2010): Audio source separation[END_REF], and a fresh dataset containing ten 3-channel mixtures of four audio sources of 10 s duration, sampled at 16 kHz. For 3-channel data we used a linear microphone array. The room reverberation time (RT) for the fresh dataset was 130 ms or 380 ms. Instantaneous mixtures are also included. Tasks T1, T2 and T3 are considered.

D2 Determined convolutive mixtures under dynamic conditions

This dataset consists of two kinds of scenarios: (1) random source activity of multiple sources in multiple static locations, and (2) a source continuously moving and overlapped with a source in a fixed or random location. The former aims to simulate a meeting scenario, where multiple talkers utter from fixed locations and their activity is unknown. The latter was specifically designed to evaluate systems able to handle dynamic variations of the mixing parameters. Due to the challenging reverberation conditions, datasets with different difficulty levels were provided (i.e. varying the source-array distance and the angular direction of simultaneously active sources). In the mixtures, two speakers are simultaneously active at most. In these datasets 4-channel mixtures are provided, and participants can decide whether using all the available channels or only a subset of them. The recordings were obtained in a real room of size (6 × 5 × 4 m) with an estimated RT of 700 ms. For both the datasets the signals were recorded by a uniform linear array of four (directional) microphones with a different spacing (of about 2 cm, 8 cm, and 18 cm) and sampled at 16 kHz. T2 and T3 are considered for this dataset.

D3 Professionally produced music recordings

According to many positive requests from the community, we decided to repeat this dataset in SiSEC2011. This dataset contains stereo music signals sampled at 44.1 kHz, including those of the dataset D3 from SiSEC2010 [START_REF] Araki | The 2010 signal separation evaluation campaign (SiSEC2010): Audio source separation[END_REF].

In addition to the 20-second snips to be separated, full-length recordings are provided as well. The mixtures were created by sound engineers, and the ways of mixing and the mixing effects applied are unknown. Task T3 is imposed on this dataset.

D4 Two-channel mixtures of speech and real-world background noise

This task aims to evaluate source separation and denoising techniques in the context of speech enhancement by merging two datasets: the dataset D3 from SiSEC2010 [START_REF] Araki | The 2010 signal separation evaluation campaign (SiSEC2010): Audio source separation[END_REF] and the CHiME corpus [START_REF] Christensen | The CHiME corpus: a resource and a challenge for computational hearing in multisource environments[END_REF]. Both datasets consist of two-channel mixtures of one speech source and real-world background noise sampled at 16 kHz. In both datasets, the spatial image of the background noise was recorded in real-world environments: a subway car, a cafeteria, or a square for the former, and a British family living room for the latter. Tasks T2, T3 and T4 are evaluated for this dataset.

All datasets include both test and development data, and the CHiME corpus in D4 also includes training data. The true source signals and source positions underlying the test data were hidden to the participants, while they were provided for the development data. The true number of speech/music sources was always available.

Evaluation criteria

Tasks T2 and T3 were evaluated via the criteria in the BSS Eval toolbox termed signal to distortion ratio (SDR), source image to spatial distortion ratio (ISR), signal to interference ratio (SIR) and signal to artifacts ratio (SAR) [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF][START_REF] Vincent | The Signal Separation Evaluation Campaign[END_REF]. In addition, version 2.0 of the PEASS toolbox [START_REF] Emiya | Subjective and objective quality assessment of audio source separation[END_REF][START_REF] Vincent | Improved perceptual metrics for the evaluation of audio source separation[END_REF] was used to assess the perceptual quality of the estimated signals for stereo data according to four performance measures akin to SDR, ISR, SIR and SAR: overall perceptual score (OPS), target-related perceptual score (TPS), interference-related perceptual score (IPS) and artifact-related perceptual score (APS).

Task T4 was evaluated by the absolute difference between the true and estimated DOAs.

Results

Despite the challenging specifications of each dataset, a remarkable participation was obtained. A total of 32 submissions were received from 18 different research centers. Many participants were involved in SiSEC for the first time, revealing a positive enlargement of the community. Tables 1 to 5 summarize the average performance obtained over the submitted algorithms. The algorithm details and all the results are available at http://sisec2011.wiki.irisa.fr/tiki-index.php. It should be noted that the presented values are the absolute values, not the improvements from the values for mixtures.

By comparison with the previous SiSEC, an unexpected high participation was observed for dataset D3. This trend seems to be in line with the recent increasing interest in NMF-based techniques, which have shown to marry well with the task of music recordings separation. The traditional dataset D1 has attracted a satisfactory amount of new participants, although the performance improvement seems to be still limited by the amount of reverberation. The datasets D2 and D4, aimed to simulate more realistic real-world scenarios, have attracted a sufficient but yet limited number of participants, probably due to the intrinsic difficulty of the data. Furthermore, the proposed algorithms do not seem to be equivalently effective in all the scenarios, which reveals that the acoustic source separation is still an open problem for real-world applications.

Note that a close analysis of each table is beyond the scope of this paper and a more detailed investigations will be discussed at the LVA/ICA 2012 conference.

Table 1 .

 1 Average performance for task T2 or T3 for instantaneous dataset D1. 2 mic: average over test & test2 datasets, 3 mic: average over test3 dataset.

	System	2 mic, 3 speech SDR ISR SIR SAR SDR ISR SIR SAR SDR ISR SIR SAR SDR ISR SIR SAR 2 mic, 3 music 2 mic, 4 speech 3 mic, 4 speech
		OPS TPS IPS APS OPS TPS IPS APS OPS TPS IPS APS OPS TPS IPS APS
	S1 [10] 13.4 25.7 21.2 14.5 16.6 27.0 23.1 20.5 8.9 17.2 15.4 9.7	-	-	-	-
		43.9 55.4 61.0 58.6 52.3 58.9 66.6 55.5 42.4 65.1 62.2 47.0	-	-	-	-
	S2 [11] 7.9	-13.6 9.7 6.9	-12.2 10.2 3.0	-	8.0 5.8 11.7 -19.1 12.6
		43.2 -61.7 25.6 40.0 -62.7 10.6 29.8 -46.7 10.7 39.7 -64.0 37.7
	S3 2	8.8	-19.8 9.4 5.9	-13.9 8.5 5.8	-16.4 6.7 8.0	-20.5 8.4
		38.5 -75.3 10.4 35.7 -68.9 16.2 35.7 -65.7 12.1 38.6 -75.5 9.2
	O1	10.8 20.1 21.7 11.1 10.4 18.0 18.8 12.5 9.1 17.6 20.0 9.3	-	-	-	-
		38.9 61.8 70.5 37.7 33.3 48.5 64.8 34.2 27.1 57.7 71.8 21.9	-	-	-	-
	O2	8.5 15.7 17.4 9.1 9.0 14.1 18.1 11.3 7.5 13.7 16.4 8.1	-	-	-	-
		24.0 29.8 72.4 20.0 30.4 28.3 69.5 21.6 22.0 20.9 70.8 13.1	-	-	-	-

Table 2 .

 2 Average performance for task T3 for convolutive dataset D1. 2 mic: average over test & test2 datasets, 3 mic: average over test3 dataset. The values are averaged over all the reverberation time.

http://mlsp2005.conwiz.dk/index.php@id=30.html

The system details can be found at the SiSEC2011 wiki.

Figure computed by averaging over an incomplete set of mixtures.

Acknowledgments:

We thank all the participants and data providers. We also thank to all the researchers who gave their opinions on the wiki and the mailinglist. Special thanks go to Mr. Antoine Liutkus, who proposed to repeat the dataset D3, and encouraged researches of the field to join in SiSEC2011.

Conclusion

This paper presented the specifications of SiSEC2011 and summarized the performance obtained over all the submissions. This time, in accordance with discussions at previous SiSECs, we carefully selected the datasets and tasks in a collaborative fashion. Ultimately, four datasets and tasks were provided which attracted many submissions from 18 research institutions. Despite some open challenges which still do not allow us to provide an unambiguous evaluation of all the submissions, we hope that SiSEC2011 will continue to represent a common platform for sharing new ideas and perspectives in the source separation research field. We believe SiSEC2011 data could be of high potential interest for many audio applications and encourage the community to use it as a reference for future evaluations.

Following the experience maturated till this campaign, new criteria seem needed for better evaluating more realistic scenarios, such as source separation Table 3. Average performance for dataset D2, "random source activity of multiple sources in multiple static locations" (top) and "a continuously moving active source overlapped with a source in a fixed or random location" (bottom). All the signals are evaluated as source signal and spatial source signal estimates. For more details see http://www.irisa.fr/metiss/SiSEC11/dynamic/main.html.

System

Source signal estimation Spatial image estimation SDR SIR SAR OPS IPS APS SDRi SIRi SARi ISRi OPSi TPSi IPSi APSi S1 [START_REF] Malek | Semi-blind source separation based on ICA and overlapped speech detection[END_REF] 3.5 9.2 7.0 30.5 69.5 11.9 2.0 6.0 7.5 3.0 29. [START_REF] Wang | On ideal binary mask as the computational goal of auditory scene analysis[END_REF] involving dereverberation or tracking of time-varying mixing conditions. Furthermore, it would be worthwhile to investigate on new objective evaluation criteria more related to the separation filter accuracy rather than to the quality of the signals itself, with the hope of minimizing the presence of outliers. With this regard, we invite all willing participants to join a continuous collaborative discussion on the future of source separation evaluation.