Joan Massich
email: joan.massich@u-bourgogne.fr

Guillaume Lemaître

Mojdeh Rastrgoo

Anke Meyer-Baese

Joan Martí

Fabrice Mériaudeau

An optimization approach to segment breast lesions in ultra-sound images using clinically validated visual cues

Keywords: Breast Ultra-Sound, BI-RADS lexicon, Optimization based Segmentation, Machine-Learning based Segmentation, Graph-Cuts

As long as breast cancer remains the leading cause of cancer deaths among female population world wide, developing tools to assist radiologists during the diagnosis process is necessary. However, most of the technologies developed in the imaging laboratories are rarely integrated in this assessing process, as they are based on information cues differing from those used by clinicians. In order to grant Computer Aided Diagnosis (CAD) systems with these information cues when performing non-aided diagnosis, better segmentation strategies are needed to automatically produce accurate delineations of the breast structures. This paper proposes a highly modular and flexible framework for segmenting breast tissues and lesions present in Breast Ultra-Sound (BUS) images. This framework relies on an optimization strategy and high-level descriptors designed analogously to the visual cues used by radiologists. The methodology is comprehensively compared to other sixteen published methodologies developed for segmenting lesions in BUS images. The proposed methodology achieves similar results than reported in the state-of-the-art.

Introduction

Breast cancer is the second most common cancer. In terms of mortality, breast cancer is the fifth most common cause of cancer death. However, it is ranked as the leading cause of cancer deaths among females in both western and economically developing countries [START_REF] Jemal | Global cancer statistics[END_REF].

Medical imaging contributes to its early detection through screening programs, non-invasive diagnosis, follow-up, and similar procedures. Despite Breast Ultra-Sound (BUS) imaging not being the imaging modality of reference for breast cancer screening [START_REF] Smith | American cancer society guidelines for breast cancer screening: update 2003[END_REF], Ultra-Sound (US) imaging has more discriminative power when compared with other image modalities to visually differentiate benign from malignant solid lesions [START_REF] Stavros | Solid breast nodules: Use of sonography to distinguish between benign and malignant lesions[END_REF]. In this manner, US screening is estimated to be able to reduce between 65 ∼ 85% of unnecessary biopsies, in favour of a less traumatic short-term screening follow-up using BUS images. As the standard for assessing these BUS images, the American College of Radiology (ACR) proposes the Breast Imaging-Reporting and Data System (BI-RADS) lexicon for BUS images [START_REF] Mendelson | Breast Imaging Reporting and Data System: ACR BIRADS -Breast Imaging Atlas[END_REF]. This US BI-RADS lexicon is a set of standard markers that characterizes the lesions encoding the visual cues found in BUS images and facilitates their analysis. Further details regarding the US BI-RADS lexicon descriptors proposed by the ACR, can be found in Sect. 3, where visual cues of BUS images and breast structures are discussed to define feature descriptors.

The incorporation of US in screening policies and the emergence of clinical standards to assess image like the US BI-RADS lexicon, encourage the development of Computer Aided Diagnosis (CAD) systems using US to be applied to breast cancer diagnosis. However, this clinical assessment using lexicon is not directly applicable to CAD systems. Shortcomings like the location and explicit delineation of the lesions need to be addressed, since those tasks are intrinsically carried out by the radiologists during their visual assessment of the images to infer the lexicon representation of the lesions. Therefore, developing accurate segmentation methodologies for breast lesions and structures is crucial to take advantage of this already validated clinical tools.

Description of the segmentation methodology

Optimization methodologies offer a standardized manner to approach segmentation by minimizing an application-driven cost function [START_REF] Cremers | A review of statistical approaches to level set segmentation: integrating color, texture, motion and shape[END_REF]. Figure 1 illustrates a generic representation of the segmentation strategy, concrete examples of its terms, applied to BUS, can be found in section 3. The overall segmentation can be seen as a three-steps strategy: (1) a mapping of the image into a discrete set of elements S, (2) the optimization stage which is formulated as a metric labelling problem, and (3) a re-mapping the labels obtained from the previous stage to produce the final delineation.

To formulate segmentation as a metric labelling problem, the image is conceived as a discrete set of elements S that need to be labelled using a label l from the labelling set L. Let W be all the possible labelling configurations of the set S, given L. Let U (•) be a cost function encoding the goodness of the labelling configuration ω ∈ W based on the appearance of the elements in S, their inner relation and some designing constraints. Then, the desired segmentation ω corresponds to the labelling configuration that minimizes this cost function, as described in Eq. ( 1).

ω = arg min ω U (ω) (1)
This goodness measure U (•) must be defined to take into account the appearance of the target region, its relation with other regions, and other designing constraints. Equation ( 2) describes this cost function as the combination of two independent costs that need to be simultaneously minimized as a whole.

U (ω) = ∑ s∈S D s (ω s ) + ∑ s∈S ∑ r∈Ns V s,r (ω s , ω r ) (2)
Where, the left-hand side of the expression integrates the so-called data term, while the right-hand side integrates the pairwise term, which is also referred to as the smoothing term. Both terms are shaped by S and evaluated in the labelling space W. In our quest to optimize the cost function U (•), it is required to define a representation for the set S, a data term D(•), a pairwise term V (•), and a proper minimization methodology.

The set S can be, in general, any discrete set representing the image (i.e. pixels, overlapping or non overlapping windows, super-pixels, etc.).

The data term D(•), given a label configuration ω ∈ W, penalizes the labelling of a particular image element or site (ω s = l) based on the data associated to s. In this manner, D s (ω s = l ) << D s (ω s = l ). Figure 2b illustrates the data cost associated to some arbitrary labelling configurations to clarify the desired effect (or behaviour) of this data term. Designing an obscure heuristic to comply with the desired behaviour of D(•) out of the box, is rather a complicated task. Therefore, an easier and cleaner approach is to design this data term D(•) with the help of Machine Learning (ML) because it provides a systematic process that is flexible enough to encode any desired behaviour based on a training stage. This concept is in fact depicted in the upper row in Fig. 1. For each site s ∈ S, features describing s are designed. Then, different optional steps can be applied to this set of features: (i) features normalization, (ii) features selection or (iii) features extraction. Finally, the data term D(•) is encoded based on ML classifiers, the features and a training step. Thus, the data term D(•) can be seen as a distance or goodness measure reflecting the likelihood for s to belong to class l.

The pairwise term V (•, •) represents the cost associated to ω s taking into account the labels of its neighboring sites, ω r , r ∈ N s . This term is usualy modeled using Markov Random Fields (MRFs) or Conditional Random Fields (CRFs). The typical form of this term, given in Eq. ( 3), is called homogenization which acts as a regularization factor favouring configurations that have a coherent labelling.

V s,r (ω s , ω r ) = { β, if ω s = ω r 0, otherwise (3)
Figure 2c shows a visual interpretation of this cost. The more fragmented is the segmentation ω, the higher is the overall pairwise term, since every boundary brings a penalization β to the total cost U (ω). In this manner, the regularization term can be seen as a post-processing or denoising stage as some sites will flip their labelling if the cost of fragmenting the regions is larger than the cost of adopting their neighbour's label.

The minimization strategy is determined by the nature of U (•) and W, since not all the minimization strategies are applicable or adequate to find ω. The size of the labelling space |W| = |L| |S| , discontinuities in U (ω) along W or the problem of local minima, additionally all the particular of all the different minimization. Need to be taken into account while choosing the most desirable minimization strategy.

BUS images segmentation using optimization

In this section, the problem of delineating structures in BUS images is defined as an optimization problem that can be solved by applying the framework presented in Sect. 2. The segmentation here proposed aims at tying a label l ∈ L (i.e. {lesion, lesion} or {chest wall, lungs, . . . , lesion}) to each element of S by simultaneously optimizing the data and pairwise terms as illustrated in Fig. 2. Choices made regarding different elements: the representation S, the data term D(•), the pairwise term V (.), and the optimizer choice are summarized in Table 1 and justified thereafter (see Fig. 1

for reference).

S is considered the result from an over-segmentation of the image using Quick-shift super-pixel [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]. The structures of the breast and their rendering when using a hand-held 2D US probe are sketched in Fig. 3a. Figure 3b illustrates the lexicon proposed by the ACR [START_REF] Mendelson | Breast Imaging Reporting and Data System: ACR BIRADS -Breast Imaging Atlas[END_REF] and used by clinicians to perform their diagnosis. Thus, our aim is to generate a set of computer vision features which is able to encode the characteristic described in the lexicon. The selected features are as follows: The relationships between the lexicon and the descriptors described previously are depicted in Table 1. More precisely, we highlight the corresponding elements of the lexicon which is encoded by each feature. A choice regarding the encoding of the data term D(•) has to be made by using a ML classifier. Support Vector Machines (SVM) classifier with Radial Basis Function (RBF) kernel is selected to determine the data model during the training stage. The pairwise term in our framework was defined as in Eq. ( 3). The optimization method used as a solver to minimize our cost function U (•) is Graph-Cuts (GC). GC, where appropriately applied, allows to rapidly find a strong local minima guaranteeing that no other minimum with lower energies can be found [START_REF] Delong | Fast approximate energy minimization with label costs[END_REF]. GC is applicable if, and only if, the pairwise term favours coherent labelling configurations and penalizes labelling configurations where neighbours' labels differ such as in our case, given by Eq. (3).

Method evaluation and comparison

The proposed methodology is evaluated using a dataset of 16 BUS images presenting a single lesion of variable extension. The size of the lesions ranges from under 1/100 to over 1/5 of the image size. The dataset is composed of cysts, Fibro-Adenomas (FAs), Ductal Inflating Carcinomas (DICs) and Inflating Lobular Carcinomas (ILCs). Every image has accompanying multi-label GT delineating all the depicted structures. This dataset is now publicly available at http://visor.udg.edu/dataset/#breast Leave-one-Patient-Out AOV:

.623 FPR: .4 FNR: .008 testing Fig. 5: Quantitative results compilation and comparison

The lack of publicly available data and source code, limits the comparison between the different methods. For this study, the results published by the other authors have been collected and expressed in terms of AOV in order to share a common metric. Further details can be found in [START_REF] Massich | Sift texture description for understanding breast ultrasound images[END_REF] and summarized in Fig. 5.

Figure 5 is divided into three main parts: (i) a table on the top summarizes the core stages of each study framework, (ii) a legend box on the right side informs our testing setup and, (iii) a comparison of the different metrics in a radial manner. An extra element is also represented in this radial representation: a blue swatch delimited by two blue dashed lines. The boundaries of this swatch correspond to the performance of some expert radiologists based on an interand intra-observer experiments carried out by Pons et al. [START_REF] Pons | Evaluating lesion segmentation in breast ultrasound images related to lesion typology[END_REF]. It is interesting to note that some methodologies outperform this swatch. A publicly available dataset should allow a better comparison in that regard.

The results point out the inherent capabilities of ML to cope with data scalability and variability, induce its usage in conjunction with larger datasets. Whereas, Active Contour Model (ACM) methodologies show its effectiveness to model the boundary in a natural manner.

For our proposed framework, the performance in terms of AOV lies within the state-of-the-art despite its final delineation limited by the capacity of the superpixels to snap the desired boundary. Figure 4 shows some qualitative results where there are limitations of labeling super-pixels when compared with handdrawn GT. Figure 4 also illustrates the influence of the pair-wise term.

Conclusions

This work presents a segmentation strategy to delineate lesions in BUS images using an optimization framework that takes advantage of all the facilities available when using ML techniques. Despite the limitation that the final segmentation is subject to the super-pixels' boundaries, the AOV results reported here are similar to those reported by other methodologies in the literature. A higher AOV result can be achieved by refining the delineation resulting from our proposed framework by post-processing it with an ACM. In this manner, the contour constraints could be applied to achieve a more natural delineation.

Fig. 1 :

 1 Fig. 1: Conceptual block representation of the segmentation methodology.

  Fig. 2: a b d e c Methodology Highlights. (a) BUS image example. (b) Superpixels' representation coloured using dataset's accompanying multi-label GT. (c) GT color code. (d) Data term: cost of labelling all sites as pectoral, lungs, adipose tissue or lesion. For illustration purposes, highly saturated colour indicates a low data cost -i.e., high confidence to assign the label associated with the color. (e) Pairwise term: labelling configurations with more boundaries produce higher pairwise term cost.

Fig. 3 :

 3 Fig. 3: Visual reference: (a) breast structures, (b) US BI-RADS lexicon

Fig. 4 :

 4 Fig. 4: Qualitative results. (a) Example 1: orignal image, super-pixels' delineations and GT. (b) Differences between GT and the delination resulting from super-pixels' boundary. (c) Ex. 2. (d) weak V (•, •) (e) strong V (•, •)

Table 1 :

 1 Design choices summary

	S	Quick-Shift super-pixels
		Background Echotexture: encoded in Appearance and SIFT-BoW
	D(•)	Echo Pattern: encoded in Appearance, Atlas and Brightness
		Acoustic Posterior: encoded in Atlas and Brightness
	V (•, •)	homogeneity as Eq. (3)
	arg min U (•)	Graph-Cuts

This work was partially supported by the Regional Council of Burgundy FEDER grant 2013-9201AAO049S02890 and by the Spanish Goverment MEC grant nb.TIN2012-3171-C02-01


Joan Massich
email: joan.massich@u-bourgogne.fr

Guillaume Lemaître

Joan Martí

Fabrice Mériaudeau

Breast Ultra-Sound image segmentation: an optimization approach based on super-pixels and high-level descriptors

Keywords: Breast Ultra-Sound, BI-RADS lexicon, Optimization based Segmentation, Machine-Learning based Segmentation, Graph-Cuts

published or not. The documents may come

Breast Ultra-Sound image segmentation: an optimization approach based on super-pixels and high-level descriptors

INTRODUCTION

Breast cancer is the second most common cancer (1.4 million cases per year, 10.9% of diagnosed cancers) after lung cancer, followed by colorectal, stomach, prostate and liver cancers. In terms of mortality, breast cancer is the fifth most common cause of cancer death. However, it is ranked as the leading cause of cancer death among females in both western countries and economically developing countries. [START_REF] Mériaudeau | Global cancer statistics[END_REF] Medical imaging plays an important role in breast cancer mortality reduction, contributing to its early detection through screening for diagnosis, image-guided biopsy, treatment follow-up and suchlike procedures. [START_REF] Smith | American cancer society guidelines for breast cancer screening: update 2003[END_REF] Although Digital Mammography (DM) remains the reference imaging modality for breast cancer screening, Ultra-Sound (US) imaging has proven to be a successful adjunct image modality. [START_REF] Smith | American cancer society guidelines for breast cancer screening: update 2003[END_REF] The main advantage of US imaging, opposed to other image modalities, lies in the discriminative power that US offers to visually differentiate benign from malignant solid lesions. [START_REF] Stavros | Solid breast nodules: Use of sonography to distinguish between benign and malignant lesions[END_REF] In this manner, US screening contributes to reduce the amount of unnecessary biopsies, [START_REF] Ciatto | The contribution of ultrasonography to the differential diagnosis of breast cancer[END_REF] which is estimated to be between 65 ∼ 85% of the prescribed biopsies, [START_REF] Yuan | Multimodality computer-aided breast cancer diagnosis with ffdm and dce-mri[END_REF] in favour of a less traumatic short-term screening follow-up using Breast Ultra-Sound (BUS) images. As the standard for assessing this BUS images, the American College of Radiology (ACR) proposes the Breast Imaging-Reporting and Data System (BI-RADS) lexicon for BUS images. [START_REF] Mendelson | Breast Imaging Reporting and Data System: ACR BIRADS -Breast Imaging Atlas[END_REF] This US BI-RADS lexicon is a set of standard markers that characterizes the lesions encoding the visual cues found in BUS images and facilitates their analysis. Further details regarding the US BI-RADS lexicon descriptors proposed by the ACR, can be found in this document in section 2.1.1, where visual cues of BUS images and breast structures are discussed to define feature descriptors.

All these facts show the interest in the medical community for incorporating US screening as part of the standard procedure in breast screening programs, [START_REF] Mendelson | Breast Imaging Reporting and Data System: ACR BIRADS -Breast Imaging Atlas[END_REF] which encourages the development of Computer Aided Diagnosis (CAD) systems using US to be applied to breast cancer diagnosis. Building CAD systems based on the clinical tools already in use (i.e. US BI-RADS lexicon) is not straight forward. Shortcomings like the location and explicit delineation of the lesions need to be addressed, since those tasks are intrinsically carried out by the radiologists during their visual assessment of the images to infer the lexicon representation of the lesions. Therefore, developing accurate segmentation methodologies for breast lesions and structures are crucial in order to develop CAD systems that can take advantage of the already validated clinical tools. This article proposes a highly modular and flexible framework for segmenting lesions and tissues present in BUS images. The proposal takes advantage of an energy-based strategy to perform segmentations based on discrete optimizations using super-pixels and a set of novel features analogous to the elements encoded by the US BI-RADS lexicon. [START_REF] Mendelson | Breast Imaging Reporting and Data System: ACR BIRADS -Breast Imaging Atlas[END_REF]

DESCRIPTION OF THE SEGMENTATION METHODOLOGY

Optimization methodologies offer a standardized manner to approach segmentation by minimizing an applicationdriven cost function. [START_REF] Cremers | A review of statistical approaches to level set segmentation: integrating color, texture, motion and shape[END_REF] Figure 1 illustrates a generic representation of the segmentation strategy here adopted to delineate breast tissues or lesions in US images. The overall segmentation can be seen as a three-steps strategy:

(1) a mapping of the image into a discrete set of elements S, (2) the optimization stage which is formulated as a metric labelling problem, and (3) a re-mapping the labels obtained from the previous stage to produce the final delineation. In order to formulate the segmentation like a metric labelling problem, the image is conceived as a discrete set of elements S that need to be labelled using a label l from the labelling set L (i.e. L = {lesion, lesion} or L = {lungs, fat, • • • , lesion}). Let W be all the possible labelling configurations of the set S, given L. Let U (•) be a cost function encoding the goodness of the labelling configuration ω ∈ W based on the appearance of the elements in S, their inner relation and some designing constraints. Then, the desired segmentation ω corresponds to the labelling configuration that minimizes this cost function, as described in Equation (1).

ω = arg min ω U (ω) (1)
It is worth to mention here that not all the minimization strategies are applicable or adequate to find ω. The convenience of a particular minimization strategy is determined by the nature of U (•) and W (see section 2.3). This goodness measure U (•) must be defined to take into account the appearance of the target region, its relation with other regions and other designing constraints. Equation (2) describes this cost function as the combination of two independent costs that need to be simultaneously minimized as a whole.

U (ω) = s∈s D s (ω s ) + s r∈Ns V s,r (ω s , ω r ) (2)
Where, the left hand side of the expression integrates the so-called data term, while the right hand side integrates the pairwise term, which is also referred as the smoothing term. Both terms are shaped by S and In general, S can be any discrete set representing the image (i.e. pixels, overlapping or non overlapping windows, etc.). For this work S is chosen to be a super-pixels representation of the image. Super-pixels can be seen as the output of an over-segmentation process or as a set of pixel collections that are contiguous and coherent with respect to a given metric. Either way, super-pixels are no overlapped irregular groups of similar connected pixels. [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF] Figure 2a illustrates one such representation S, applied to a BUS image example. The superpixels are coloured according to the image's Ground Truth (GT). For the rest of this work, S is considered to be the super-pixels resulting from an over-segmentation of the image using Quick-shift. [START_REF] Massich | Deformable object segmentation in ultra-sound images[END_REF] Bear in mind that given an unseen BUS image, the ultimate goal is to represent the image as a set of superpixels and infer the appropriated labelling for each of them. This goal requires to define: a data term, a pairwise term, and a proper minimization methodology.

The data term

Given a label configuration ω ∈ W, the data term penalizes the labelling of a particular image element or site (ω s = l) based on the data associated to s. In this manner, D s (ω s = l ) << D s (ω s = l ). Figure 2b illustrates the data cost associated to some arbitrary labelling configurations to clarify the desired effect (or behaviour) of this data term (fig. 2a shows the GT of each site s). Notice that the labelling configurations ω used in fig. 2b have the particularity that all sites share the same label, ω ∈ {ω s = l, ∀s ∈ S}.

Designing an obscure heuristic to comply with the desired behaviour of D(•) out of the box, is rather a complicated task. Therefore, an easier and cleaner approach is to take advantage of Machine Learning (ML) techniques to design this data cost in a systematic manner based on a training stage. The idea is to generate a data model for each label (or class) in L from training samples, and let D(•) be a distance or goodness measure reflecting the likelihood for s to belong to class l. ML offers a systematic but flexible approach to customize D(•). Different features to represent the data, custom construction of the data models by using different classifiers, training techniques or including arbitrary constraints; can be used to achieve the desired data term without changing the overall scheme.

Despite details regarding the construction D(•) are out of the scope of this report, the rest of this section 2.1 summarizes this process. For further details the reader is referred to Massich et al. [START_REF] Massich | Deformable object segmentation in ultra-sound images[END_REF] The usage of ML as part of the proposed framework to determine D(•) is represented at the upper side of the diagram in fig. 1

BUS features to build the data term

Figure 3 is a three-parts illustration to graphically summarize the visual cues that can be found in BUS images and their incorporation to the data term. The structures of the breast and their rendering when using a handheld 2D US probe are shown in fig. 3a. Figure 3b illustrates the lexicon proposed by the ACR. [START_REF] Mendelson | Breast Imaging Reporting and Data System: ACR BIRADS -Breast Imaging Atlas[END_REF] Whereas fig. 3c relates the visual cues to the following features:

Appearance Based on the multi-labelled GT, a Median Absolute Deviation (MAD) histogram model for every tissue label is built. The Appearance feature is computed as the χ 2 distance between histogram of s and the models.

Atlas Based on the multi-labelled GT an atlas is build to encode the label likelihood based on the location of s.

Brightness Takes an intensity descriptor of s (i.e: mean, median, mode) and compares it with some intensity markers of the set S such as the minimum intensity value, the maximum, its mean, etc.

Self-Invariant Feature Transform (SIFT)-Back-of-Features (BoF) s is represented as the occurrences of a SIFT dictionary of 36 words. [START_REF] Massich | Sift texture description for understanding breast ultrasound images[END_REF] In order to incorporate multi-resolution, each super-pixel is group with its adjacent super-pixels such that s = {s ∪ N s }, the features are recalculated using s and concatenated to the original feature descriptor of s. This operation can be repeated several times.

Notice that not all the visual cues compiled by the lexicon are suitable to build the data term. Mass shape, orientation and margin cues do not characterize an individual super-pixel but a group of contiguous super-pixels sharing the same label. Characterising super-pixels in those terms imply that super-pixels should be image elements of the same order as the lesions. This is not the case, since the bottom-line is to aggregate small regions to form the segmentation. The lesion boundary cue falls at the other extreme. To encode this visual cue as data term feature implies to use smaller super-pixels in order to be fully contained within the halo, limiting the discriminatory power of other features that need larger regions to build higher-level descriptors.

In this manner the B.Echo-texture, is mainly encoded by the SIFT-BoF but also through the Appearance feature since elements would have the same model. Acoustic Posterior is mainly captured by the Brightness feature, but the Atlas also brings crucial information to compensate for intensity inhomogeneities caused by signal attenuation present in the most posterior parts of the image. The Echo Pattern cue, encodes the echogenity of a region with respect to the adipose tissue. This can be found at the anterior part of the image with its echogenity close to the middle of the spectrum. Appearance, Brightness and the Atlas account for such information.

The data term construction

For the work here reported, a Radial Basis Function (RBF)-Support Vector Machines (SVM) classifier is used to determine the data model during the training stage. During the testing stage, D s (ω s = l) corresponds to the distance between the testing samples and the model associated to l as the SVM classification reward.

The pairwise (or smoothing) term

The pairwise term represents the cost associated to ω s taking into account the labels of its neighbour sites, ω r , r ∈ N s . This term models a Markov Random Fields (MRFs) or a Conditional Random Fields (CRFs). The typical form of this term, given in eq. ( 3), is called homogenization which acts as a regularization factor favouring configurations that have a coherent labelling.

V s,r (ω s , ω r ) = β, if ω s = ω r 0, otherwise (3)
Figure 2c offers a visual interpretation of this cost. The more fragmented is the segmentation ω, the higher the overall pairwise term; since every boundary brings a penalization β to the total cost U (ω). In this manner the regularization term can be seen as a post-processing or denoising stage since that some sites will flip their labelling if the cost of fragmenting the regions is larger than the cost of adopting their neighbour's label.

More sophisticated smoothing terms, where boundaries have different penalization based on site relations in S and image information as illustrated in fig. 2c, are also naturally handled by the proposed framework. [START_REF] Massich | Deformable object segmentation in ultra-sound images[END_REF]

Searching the best labelling configuration

Once defined U (ω) so that the cost for a particular labelling configuration ω can be computed, the problem of finding ω corresponding to the global minimum of the space W of all possible labelling configurations needs to be faced. This problem falls into the category of NP-hard problems. Moreover, due to limitations in building U (•) such as noise, training policies, etc. there are no guarantees that the global minimum ω corresponds to the true labelling.

Nevertheless, there is a large body of literature proposing methodologies to find suboptimal solutions to the problem trading-off between time of convergence and accuracy of the solution reached. Szeliski et al. [START_REF] Szeliski | A comparative study of energy minimization methods for markov random fields with smoothness-based priors[END_REF] conducted an exhaustive review in terms of solution quality and runtime of the most common energy minimization algorithms used in Computer Vision (CV), such as Iterated Conditional Modes (ICM), Simulate Anealing (SA) or Graph-Cuts (GC).

The minimization strategy used for this work is GC. This technique was initially introduced to solve CV applications by Boykov et al. [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] Soon after its introduction, it becomes the minimization technique of choice for CV problems. Since, when GC is applicable, it allows to rapidly find a strong local minima guaranteeing that no other minimum with lower energies can be found. [START_REF] Delong | Fast approximate energy minimization with label costs[END_REF] GC is applicable if, and only if, the pairwise term favours coherent labelling configurations and penalizes labelling configurations where neighbours labels differs; such is our case, given in eq. (3).

METHOD EVALUATION AND COMPARISON

A 16 BUS images dataset with accompanying multi-label GT delineating all the structures present in the images has been used to evaluate the proposed methodology for lesion segmentation application. Every image in the dataset presents a single lesion with variable extension. The size of the lesions ranges from under 1/100 to over 1/5 of the image. The dataset composed of cysts, Fibro-Adenomas (FAs), Ductal Inflating Carcinomas (DICs) and Inflating Lobular Carcinomas (ILCs).

Figure 4 shows qualitative results, results whereas the quantitative results from the best configuration are reported as a table in fig. 5b. [START_REF] Massich | Deformable object segmentation in ultra-sound images[END_REF] Notice that for the proposed framework, the performance in terms of Area Overlap (AOV) is limited by the capacity of the super-pixels to snap the desired boundary. Figure 4b shows how the delineation resulting from a proper labeling of the super-pixels differs from the GT. Despite the fact that a False Positive Ratio (FPR) of 0.4 seems significant, based on our experiments most of the images produce no False Positive (FP) lesions. However, images with FP lesions, are likely to produce more several of them as shown in fig. 4c-e. The amount of FP lesions can be trimmed by applying a higher cost in the pairwise therm (compare fig. 4d and fig. 4e. Nevertheless, a severe increasing of the pairwise cost also increases the False Negative Ratio (FNR) since some lesions are missed due to over-smoothing. [START_REF] Massich | Deformable object segmentation in ultra-sound images[END_REF] The situation of having a larger FNR is less desirable than reducing the FPR. The FNR reported in fig. 5b is caused by an image within the dataset that its lesion is fully contained in a single super-pixel and still around 20% of this super-pixel's area is healty tissue. Due to the lack of publicly available data (and code) there is no manner to perform a methodology comparison further than compiling the results reported in the literature. The table in Figure 5a complies the evaluation reported by the authors of the most relevant methodologies found in the literature. [START_REF] Massich | Deformable object segmentation in ultra-sound images[END_REF] Details about the methodologies proposed in the literature can also be found in the aforesaid table. Specifically it is detailed the category of technique been used for detecting the lesions, segmenting it and post-process the delineations (if any). The studied categories are: ML, Active Contour Model (ACM) and others. The iconography used in fig. 5a also illustrates if those stages are treated as independent and connected in a daisy-chain fashion, or otherwise the stages are addressed in an atomic manner.

Figure 5c renders the information present in fig. 5a andfig. 5b in a visual manner compare all the results at once. The methodologies arranged in a radial fashion and grouped by its most representative technology category. In red it can be found a small, medium and large categorization of the dataset reported for testing. The concentric circles represent AOV. The blue line correspond to the AOV results reported in the literature, whereas the black line indicates the AOV our framework scored in our testing. An extra element is also represented in fig. 5c as blue swatch delimited by two blue dashed lines. The boundaries of this swatch correspond to performance of expert radiologists in terms of AOV based on an inter-and intra-observer experiment carried out by Pons et al. 14 * .

When comparing the results it is clear the inconvenience of unexciting public data, since several of the results outperform the manual delineations studied in. [START_REF] Pons | Evaluating lesion segmentation in breast ultrasound images related to lesion typology[END_REF] It can also be seen that the category tested in larger datasets is ML, whereas ACM lead to better segmentations since the lesion boundary is easier to model in ACM compared to ML based techniques.

CONCLUSIONS

This work presents a segmentation strategy to delineate lesions in BUS images using an optimization framework that takes advantage of all the facilities available when using ML techniques. Despite the limitation that the final segmentation is subject to the super-pixels' boundaries, the AOV results here reported are similar to those reported by other methodologies in the literature. A higher AOV result can be achieved by deforming the delineation resulting from the proposed framework using this to initialize a second post-processing step based on ACM. In this manner the contour constrains could be applied to achieve a more natural delineation.

Figure 1 :

 1 Figure 1: Conceptual block representation of the segmentation methodology

Figure 2 :

 2 Figure 2: Methodology terms interpretation

  , which can be divided into two blocks: (a) the features to represent the samples, and (b) the tools to encode D(•) based on the features and the training.

Figure 3 :

 3 Figure 3: Visual reference of breast structures and visual cues used for standard BUS image assessment and diagnosis.

Figure 4 :

 4 Figure 4: Qualitative results. (a) Example 1: orignal image, super-pixels' delineations and GT. (b) Differences between GT and the delination resulting from super-pixels' boundary. (c) Ex. 2. (d, e) Pairwise term influence to FPR, FNR

  (a) BUS images lesion segmentation strategies compiled from the bulk of the literature: reported quantitative results and methodology highlights. technology used for: detection segmentation post-processing S: Quick-Shift super-pixels arg min(U (•)): Graph-Cuts D(•):

Figure 5 :

 5 Figure 5: Quantitative results compilation and comparison

This work was partially supported by the Regional Council of Burgundy.


Mohamed Belkacemi
email: mohamed.belkacemi@u-bourgogne.fr

Christophe Stolz

Alexandre Mathieu

Guillaume Lemaître

Joan Massich

Olivier Aubreton

Non Destructive Testing based on a Scanning-From-Heating approach: Application to non-through defect detection and fiber orientation assessment

Keywords: Non Destructive Testing, Active Thermography, 3D scanning, Fiber orientation, Non-through defect detection

come

Introduction

Active infrared thermography becomes one of the essential techniques for Non-Destructive Testing (NDT). Approaches like Pulse Thermography 1 Pulse Phase Thermography 2 or Lock-in Thermography [START_REF] Wu | Lock-in thermography for nondestructive evaluation of materials[END_REF] have proved to be effective to characterize defects. In recent years, an evolution of NDT systems by adding 3D measurement can be noted. Fernandes et al. [START_REF] Fernandes | Fiber orientation assessment in complex shaped parts reinforced with carbon fiber using infrared thermography[END_REF] reconstruct carbon fiber objects using 3D scanner and determine the orientation of the fiber with a laser and an infrared camera. Subsequently, the thermal information is mapped to the reconstructed 3D model. Similarly, Oswald-Tranta et al. [START_REF] Oswald-Tranta | Fusion of geometric and thermographic data for automated defect detection[END_REF] detect crack in non-magnetic materials based on the fusion between the infrared and the 3D data, obtained through an inductive NDT system and a conventional stereo vision system, respectively. To quantify the energy losses in buildings, Vidas et al. [START_REF] Vidas | 3d thermal mapping of building interiors using an rgb-d and thermal camera[END_REF] map the thermal information on the 3D data using a Microsoft Kinect and a thermal camera. To conclude, it can be noted that most of the proposed prototypes are actually using two systems: a 3D scanner and a non-destructive inspection system in the infrared range. These hybrid systems in which distinct steps of 3D reconstruction and NDT inspection suffer from drawback such as errors during the data fusion.

Scanning-From-Heating (SFH) offers an alternative. This approach is initially used to reconstruct objects into 3D. Indeed, Eren et al. [START_REF] Eren | Scanning from heating: 3d shape estimation of transparent objects from local surface heating[END_REF] achieve impressive results using SFH in which the geometry of glass or plastic transparent objects is estimated based on active triangulation. The transparent object placed on a gliding stage is locally heated using a CO2 laser source of 10.6 µm wavelength. The emitted radiation is acquired by an Infra-Red (IR) camera and the 3D position of the heated point is recovered thanks to a prior calibration, with an accuracy in the order of 100 µm. Bajard et al. 8 have extended this method for metal objects which have different radiative and conductive properties. The CO2 laser source is replaced by an Nd :YAG Laser of 1.06 µm wavelength, since the metals absorption coefficient is 4 times greater at 1.06 µm than at 10.6 µm. The results of the specular objects digitization are shown in Fig. 1. Other examples of 3D reconstruction are presented for black plastic, composite material, and ceramic using the same system.

In this work, we present an approach in which SFH is extended to perform both 3D reconstruction and NDT. We aim at detecting non-through defect and estimate the orientation of fibers using punctual stimulation. The proposed approach significantly reduces the cost of the hardware setup and avoids data fusion errors. Furthermore, a simulation is carried out to demonstrate that the non-defective area can be localized using the thermal radiation disturbance.

The rest of the paper is organized as follows: Section 2 summarizes related works regarding NDT based on punctual stimulation. The proposed method is explained in Sect. 3 while Sect. 4 presents experimental results. Conclusion and avenues for future directions are drawn in Sect. 5.

Related Work

Generally, active thermography approaches use a uniform heating over the entire object-surface and the data are processed assuming the 1D model of the heat transfer. However, several works tackle this problem using local heating. Burrows et al. 9 detect crack defects on object surface using a laser beam. The heat distribution is disrupted by the crack which highlight cracks. In the case of non-through defects, Hammiche et al. 10 make use of a Thermal Scanning Microscope where a probe is put into contact with the surface of the sample. Thermal contrast between the excitation and the response of the material is used to build a contrast image at each scanned point of the object. Along the same lines, Ermert et al. 11 develop a non-contact detection approach based on a Scanning Electron Microscope. A focused and modulated electron beam excites the object surface and the infrared radiations are captured by a pyroelectric detector. However, the surface controlled is restricted.

Fiber orientation can be estimated using pulsed thermal ellipsometry technique 12 which consists in heating the part using a laser beam. The inspected part is heated by a laser beam and due to the anisotropy of the material, the observed thermal pattern becomes elliptical. Subsequently, The orientation of the fiber is deduced from the ellipse orientation.

3 Non-Destructive Testing based on punctual excitation

Non-through defect detection

In this section, the concept used to detect defect is presented and validated through simulation.

Then, a specific image processing algorithm is developed to analyze the thermal evolution of each heated point.

Simulation

A 3D transient analysis is developed on COMSOL c software using Finite Element Method. The test object is a steel plate with 10 mm thickness with the following properties: one face is flat while the other is degraded by non-through holes of 10 mm diameter drilled at different depths. A thermal simulation including thin insulating layer subjected to laser irradiation is carried out. A conductive heat transfer within the metal medium is simulated through the time dependent partial differential equation such as:

ρC p ∂T ∂t = ∇ [k∇T ] + Q , (1)
where ρ is material density, C p is specific heat capacity, T is temperature, t is time, k is thermal conductivity and Q is a heat source term. In our simulation, Q has been set to 0 and the metal medium is considered as opaque. The laser heat source is considered as an inward heat flux and characterized by a Gaussian distribution such as:

q 0 = 2 P L πr 2 0 exp -2 x 2 + y 2 r 2 0 , (2)
where P L is the laser power, r 0 is the laser beam waist, x and y are spatial coordinates.

Furthermore, the black paint layer is considered as a dielectric material so that the energy provided by the interaction between the laser beam and material involves an absorption coefficient α, near 1. Transient conduction analysis of heat transf er module has been selected to solve heat equation. The 2D geometry of the problem is represented in Fig. 2a where d is the defect depth, D is the defect diameter. Boundary condition imposed to all surfaces is in form of the following equation:

n(k∇T ) = αq 0 + h(T ext -T ) , (3)
where n is the normal vector, α the absorption coefficient, q 0 is the inward heat flux on top surface and is equal to 0 on other surfaces, h is a coefficient representing heat transfer losses due to convective heat flux with surrounding atmosphere. h is equal to 5 W m -2 K -1 on all surfaces, except on the surfaces of the defect where it is equal to 0 to simulate thermal insulation. T ext is the surrounding atmosphere temperature and it is assumed that the radiation heat transfer losses are neglected.

For meshing, free tetrahedral meshes are selected with choice of a custom element size. The maximum element size is equal to 0.032 mm with a maximum element growth rate equal to 1.35.

The smallest element size is obtained near to the origin (0, 0, 0) centrally located on top of the surface (see Fig. 2c). A fine meshing is constructed with a number of elements ranging from 75, 000 to 225, 000 due to the presence of a 0.1 mm layer on top of the simulated object surface. Furthermore, the system to solve consists of a number of degrees of freedom ranging from 111, 000

to 324, 000 and depends of the chosen geometry. Record time stepping is selected to 0.1 s in a range from 0 s to 2 s. Geometric Multigrid GMRES iterative linear solver is used. Backward Differentiation Formula (BDF) implicit method is used to help solver to converge.

The interaction of the laser beam with the solid surface produces a thermal wave propagating through the solid media. In the case of a media without any defects and assuming a semi-infinite and isotropic solid, the thermal response is defined by the shape of a symmetric hemisphere 13 as shown in Fig. 2b. This response is modified due to the presence of a defect close to the excited area.

The simulation results are depicted in Fig. 3a-3b-3c by representing the temperature distribution in kelvin. It can be noted that the temperature intensity increases when the laser heats the defective area. Subsequently, the maximum temperature is reached for defects nearer to the surface. Fig. 3d-3e-3f represent the thermal response using another color-map. It can be observed that the amount of the thermal response increases due to the presence of inhomogeneities in the solid media.

Defect detection algorithm

The object surface is punctually scanned using a laser beam and at each of these scans, a thermal image is acquired (see Fig. 7c). Defect localization is characterized by a change of the heat distribution and therefore can be detected with an appropriate image processing algorithm in the thermal images. Our developed framework is depicted in Fig. 4 and can be subdivided into three stages: (i)

pre-processing, (ii) segmentation, and (iii) detection.

Pre-processing Each acquired thermal image is enhanced by successively subtracting the background, reducing the noise, and stretching the contrast. The background to be subtracted consists in the acquisition of a thermal image before laser stimulation. Additionally, a median filtering is applied to the image to reduce the noise and the denoised image is enhanced through contrast stretching.

Segmentation Then, the enhanced image is segmented using the region growing segmentation algorithm. 14 The initial seed point is automatically placed on the center of the thermal response corresponding to the maximum temperature in the image. Then, the region grows if the difference of temperature values of the neighboring pixels and the region's mean is less than a given region membership criterion. The pixel with the smallest difference is assigned to the segmented region.

Detection The detection is subdivided into two stages: (i) merging and (ii) binarization. All segmented images are combined together by creating a matrix s such that the number of white pixels of a given segmented image is assigned to the location (i, j), which corresponds to the location of the laser stimulation. Subsequently, to limit the influence of the emissivity change and detect robustly the defective area, the median s is computed for each Region Of Interest (ROI) with a size (40 × 40) and compared to each value s(i, j) inside this ROI in order to binarize s:

s m (i, j) =          1 if s(i, j) ≥ K × s 0 if s(i, j) < K × s (4)
where K > 1 and s m (i, j) is the resulting binary mask of the defects detection. If K 1, the defect detection is more robust to noise. However, only shallowest defect are detected; if K → 1, the defect detection is subject to noise. However, deeper defect are detected.

Fiber orientation assessment

The fiber orientation is estimated based on the pulsed thermal ellipsometry approach. 12 Indeed, the thermal response of a polymer matrix composite material stimulated by a laser beam is elliptic due to the propagation of the heat through the fiber. The framework developed to infer the fiber orientation is depicted in Fig. 5. Thermal images are acquired with an IR camera and are pre-processed with the following steps: (i) background subtraction, (ii) a segmentation using region growing, and

(iii) an edge detection using Canny edge detector. 15 Three methods are commonly used to fit an ellipse: (i) algebraic fitting, 16 (ii) orthogonal least squares fitting, 17 and (iii) Maximum Likelihood algorithms. 18 The two latter methods are suitable for 3D data, whereas the first method can only be applied to 2D problems. Algebraic fitting, however, offers a low computational complexity and is suitable in our case. In this method, the ellipse is formulated by the general conic polynomial:

F (a, X) = aX = 0 (5)
with

a = A B C D E F X = x 2 xy y 2 x y 1
The fitting is performed by a least square minimization of the Euclidean distances such that:

arg min a n i=1 Ax 2 i + Bx i y i + Cy 2 i + Dx i + EY i + F 2 . ( 6
)
The following constraint yields to one solution which corresponds to an ellipse:

4 Experiments and results

Experimental setup

Figure 6 illustrates the complete scanning system setup. The system is composed by a FLIR 645 IR camera which captures the thermal radiation and an Nd :YAG Laser system to stimulate the object surface. The FLIR 645 infrared camera provides a sensitivity range from 7.5 to 14 µm. The Nd :YAG Laser system produces a step signal of 1.5 s period and the power of the laser beam is set up at 1.5 W. A galvanometric mirror is used to control the position of the beam, which facilitates the object scanning with a 1 mm resolution. During the experiment regarding the fiber orientation assessment, the motorized rotary stage ZABER T-RS60A is used.

Experimental temperature profile

In order to validate the simulation results presented in Sect. 3.1.1, a set of thermal images are acquired. In these images, a defective area with a defect of 0.5 mm below the surface and nondefective areas are excited with the laser beam (see Fig. 7). As previously predicted in the simulation section, the thermal response of a defective-area is different from the non-defective area, allowing non-through defects detection.

Non-through defect detection

Two type of material are used for the non-through defect detection: (i) steal plate and (ii) aluminum plate. Both type of plates have non-through artificial defects, i.e. circle, square, triangle, and cross. Defects are located between 0.5 mm and 2 mm below the observed surface for the aluminum plate and between 1 mm and 3 mm in the case of the steel plate. The images are processed with the image processing algorithm presented in Sect. 3.1.2. A 4-neighbors is considered during the region growing algorithm and a membership criterion of 0.2, while the matrix s m is generated with

K = 1.4.
The results are depicted in Fig. 8 by superimposing with a color map on the raw image (i.e., red pixels are the non-through defects).

Figure 8a shows a steel plate with circular non-through defects, in which defects have a diam-eter of 9, 7, 5, and 3 mm, from left to right. Each defect is positioned at a different depth 1, 2, 2.5, and 3 mm, from bottom to top. Only defects with depth less than or equal to 2.5 mm are detected.

However, the defects with a diameter less than 5 mm are not detected as shown in Fig. 8b. Besides, for the aluminum plate all defects are detected and the shape of the defect can be distinguished (see Fig. 8d).

To study the repeatability of the technique, our framework is evaluated on steel plate with 16 identical circular defects located at 1 mm depth and a 10 mm diameter. The result of this experiment is depicted in Fig. 9b, in which all defects are successfully detected.

Fusion of defect localization and 3D data

To assess the quality of the fusion between the defect localization and the 3D data, a steel object is scanned (see Fig. 10a) to acquire a set of thermal images. The steel object contains non-through circular defects in the planar part with a depth of 1 mm. The scanned area corresponds to the highlighted red delimitation depicted in Fig. 10a. A similar method as proposed by Bajard et al. 8 is used to obtain the 3D information and the defects are detected as in Sect. 3.1.1. The coordinate system is the same for the 3D and defects information due to the fact that these data are acquired with the same system. Therefore, the fusion between these data is straightforward, which overcome the matching problem faced by other methods as stated in the introduction section.

Fiber orientations assessment

In this section, we evaluate our framework at detecting the fiber orientation of carbon fiber composite material (see Fig. 11a). The fiber orientation is estimated using the methodology presented in Sect. 3.2. The initial experimental setup is as follows: firstly, the material is heated with the laser beam, the thermal radiation is recorded, and the orientation of the fiber is inferred. Secondly, the experiment is repeated by rotating the carbon plate at a 90 • position. The results of the ellipse fitting are depicted in Fig. 11b-11c-11d. The algebraic fitting qualitatively shows an accurate fitting as shown in Fig. 11d, in which the blue dots represents the raw data, the green curve is the fitted ellipse, and the red line represents the ellipse major axis orientation. As expected, the ellipse orientation follows the material rotation as depicted Fig. 12.

The previous experiment are repeated acquiring 4 measurements for different orientations ranging from -10 • to -60 • with a step of 10 • . The absolute orientations are estimated and the mean and standard deviation are computed as presented in Table 1. The average relative orientation is equal to 10.8 • with a standard deviation of 1.7 • and can be compared to the step which is set to

10 • .
Table 1 Fiber orientation evaluation.

Initial position -10 • -20 In this paper, a system based on the SFH approach is proposed for detecting non-through defects and fiber orientation. In the one hand, experiments conducted on steel and aluminum plates showed that non-through defects until 2.5 mm are detected. On the other hand, a test was carried out in order to investigate the influence of the depth and the size of the defects on the defect detection. In addition, accurate fiber orientation assessment is achieved with an experimental standard deviation of 1.7 • . The use of the SFH for NDT application and the 3D digitization offers the prospect to imagine an industrial system that can give a complete solution for quality control process. As future works, tests will be done on different specimen to explore the influence of material properties, defect size and depth on the detection method. Detecting defect on complex objects is a though challenge in active thermography. Our technique has proven to be effective for non-through defect detection on planar object. The next step is to investigate the possibility of detecting defects on complex geometry objects. Moreover, kinetic analysis of the thermal radiation could allow to estimate the depth of the defect.

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)

Fig 1 Fig 2

 12 Fig 1 Examples of 3D digitization obtained by SFH approach: (a), (b), (c) and (d) metallic specular surfaces 19 -(e) and (f) black plastic -(g) and (h) composite material -(i) and (j) ceramic -(k) and (l) glass transparent object. 20

Fig 3 Fig 4 Fig 5

 345 Fig 3 (a) and (d) Thermal response, observed on a non-defective area -(b) and (e) Thermal response, observed on an area with defect of 1 mm below surface -(c) and (f) Thermal response, observed on an area with defect of 0.5 mm below surface.

Fig 6

 6 Fig 6 Scanner prototype.

Fig 7

 7 Fig 7 (a) Front side of the inspected material -(b) Back side of the inspected material -(c) Form of the thermal response of the material in the area with and without defect.

Fig 8

 8 Fig 8 (a) and (c) Steel and aluminum objects -(b) and (d) Corresponding detection map.

Fig 12

 12 Fig 9 (a) Steel objects -(b) Corresponding detection map.

  • -30 • -40 • -50 • -60 •

	5 Conclusion and future works
		42.2	31.3 23.3 12.2	0.2 -10.9 -22.1
	Measurements	42.8	31.1 22.5 13.9	0.4 -10.9 -21.7
		42.8	30.3 23.1 13.1 -0.9 -10.7 -22.1
		43.7	31.8 22.9 13.1	0.4 -10.9 -22.1
	Average	42.9	31.1 23.0 13.0	0.0 -10.9 -22.0

Acknowledgments

This work is supported by the Regional Council of Burgundy. We are extremely grateful to members of laboratory LE2I for day to day help and collaboration.


Antoine Perrier

Marek Bucki

Vincent Luboz

Nicolas Vuillerme

Yohan Payan

A Perrier Abc

3D musculo-skeletal finite element analysis of the foot kinematics under muscle activation with and without ankle arthrodesis

Keywords: Foot model, Finite Element, Musculoskeletal, Motion Analysis, Arthrodesis

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

The choice between arthrodesis and arthroplasty in the context of advanced ankle arthrosis remains a highly disputed topic in the field of foot and ankle surgery. Arthrodesis, however, represents the most popular option. Biomechanical modeling has been widely used to investigate static loading of cadaveric feet [START_REF] Valderrabano | Kinematic changes after fusion and total replacement of the ankle: part 3: Talar movement[END_REF] as well as consequences of arthrodesis on bony structures. Although foot kinematics has been studied using motion analysis, this approach lacks accuracy in capturing internal joints motion due to limitations inherent to external "marker sets" and the fact that it imposed the foot to be considered as a rigid solid. The consequences of arthrodesis on kinematics of the unloaded foot are not well understood although it is of crucial importance during the swing phase and at heel contact. Investigating ankle mobility during muscle contraction with and without arthrosis could explain how the motion is produced by extrinsic muscles activations affected by an arthrodesis. This study aims at defining if a biomechanical model with Finite Elements (FE) could help arthrodesis understanding.

Methods

Modeling

A multi-body musculoskeletal model of the foot has been developed using the ArtiSynth simulation platform [START_REF] Lloyd | ArtiSynth: A Fast Interactive Biomechanical Modeling Toolkit Combining Multibody and Finite Element Simulation[END_REF]. The model features rigid bodies modeling the bones and joint interactions as well as FE modeling the soft tissues of the foot.

• Musculoskeletal model A high-resolution CT volume acquired on the unloaded right foot of a patient was used to segment the 30 bones and reconstruct their 3D shapes. Rigid body contact constraints have been implemented in the Artisynth framework to model joint interactions. All 33 joint motions are further constrained by 210 ligaments modeled by cables and inserted on the bones using CT images of the patient. The foot model does not rely on idealized joint models but solely on the modeling of contacts between bones and of ligaments attached to the bones in order to guide the foot kinematics. The Aponeurosis has been modeled using five ligaments linked by transverse structures. 1) The first simulated movement features the activation of toes and hallux extensor muscles alone, followed by a return to a resting position.

2) The second involves the contraction of the tibialis anterior muscle alone, followed by a return to a resting position.

3) The third simulated movement is driven by the synergistic contraction of the extensors and the tibialis anterior. Plantar flexion activation patterns initiate the foot movement by first contracting muscles that allow the plantar flexion of the ankle (Tab 1, bottom lines):

1) The first simulated movement features the activation of toes and hallux flexor muscles alone, followed by a return to a resting position.

2) The second involves the contraction of the tibialis posterior muscle alone, followed by a return to a resting position.

3) The third simulated movement involves the contraction of the triceps, followed by a return to a resting position. 4) The fourth simulated movement is driven by the synergistic contraction of all those ankle flexor muscles.

Evaluation

The 3D angle called hereafter "Foot motion", formed by the mass center of the tibia, the mass center of the talus and the center of the second metatarsal head, is considered as an anatomical axis of the foot. In the simulations performed with arthrodesis, only the first 3D angle, called hereafter "Ankle fusion -Foot motion", was computed given that the tibia-talus angle was set to 117° by the simulated joint fusion. Angle measurements were also performed in 2D after projecting the reference points on the three anatomical planes. Movements were split into:

• Dorsal/plantar flexion in the sagittal plane,

• Pronation/supination in the frontal plane, and • Abduction/adduction in the horizontal plane. These three categories are made to facilitate a clinical interpretation of the results.

Results and discussion

The combined action of the anterior muscles induces a 3D dorsal flexion with peak values: Foot motion = 42.2°. In this activation pattern, the arthrodesis limits dorsal flexion to Ankle fusion -Foot motion = 8.7.

The combined action of the posterior muscles induces a 3D plantar flexion with peak values: Foot motion = 32.7°. With the arthrodesis, 3D plantar flexion is Ankle fusion -Foot motion = 14.4°. By simulating a virtual arthrodesis, the biomechanical model makes it possible to explore the relationships between adjacent joints and to quantify their respective amplitudes within a foot motion. Hence the contraction of the TA muscle alone (2 nd movement) with the arthrodesis makes it possible to compute the mobility of the navicular on the talus. Indeed, 2.8° movement (Ankle fusion -Foot motion angle in the sagittal plane during dorsal flexion) is produced by the navicular alone since the talus is fixed. The missing Ankle fusion -Foot motion amplitude in the sagittal plane (i.e. 8.7° -2.8° = 5.9°) is provided through the Cuneo-metatarsal and Cuneo-navicular joints, which are mobilized by the synergistic contraction of anterior muscles Tab 1: 3D Foot motion with and without arthrodesis under muscle activation. Extensus (E), Flessor (F), Hallucis Longus (HL), Digitorum Longus (DL),Tibialis Anterior (TA), Posterior (TP)

Conclusion

A patient-specific model of the foot with and without simulated arthrodesis, driven by muscle activations reproducing the swing motion of the foot, could be used for surgical planning to find the optimal setting of the tibio talar angle in order to optimize swing phase during gait.

  Finally, 15 Hill's model muscles have been positioned according to their anatomical course and can be independently activated (Fig 1). • Finite Element model Soft tissues have been modeled by a FE mesh comprising three sub-domains representing skin, fat and muscles tissues. The resulting mesh has 78,161 elements (21,900 hexahedrons, 17,319 pyramids, 22,304 tetrahedrons, and 16,638 wedges) and 47,242 nodes. Neo-Hookean materials were used to account for large deformations. Each soft-tissue sub-domain implements the relevant mechanical properties, namely: E=200kPa, ν=0.485 for the skin; E=60kPa, ν=0.495 for muscles; and E=30kPa, ν=0.49 for fat tissues (Fig 2). All the modeled elements are subject to gravity.

Fig 1 :

 1 Fig 1: Musculoskeletal (left) and finite element (right)parts of the foot model.2.2 SimulationA simulation using hill's muscle activation patterns for dorsal and plantar flexion of the ankle with and without arthrodesis was carried out. Muscles are activated by shortening their length. The mobility of the talus bone and the foot was assessed using kinematic measurements. Dorsal flexion activation patterns initiate the foot movement by first contracting muscles that enable the dorsal flexion of the ankle (Tab 1, top lines):

Fig 2 :

 2 Fig 2: Dorsiflexion under three different muscle activations.2.3 EvaluationThe 3D angle called hereafter "Foot motion", formed by the mass center of the tibia, the mass center of the talus and the center of the second metatarsal head, is considered as an anatomical axis of the foot. In the simulations performed with arthrodesis, only the first 3D angle, called hereafter "Ankle fusion -Foot motion", was computed given that the tibia-talus angle was set to 117° by the simulated joint fusion. Angle measurements were also performed in 2D after projecting the reference points on the three anatomical planes. Movements were split into:• Dorsal/plantar flexion in the sagittal plane,• Pronation/supination in the frontal plane, and • Abduction/adduction in the horizontal plane. These three categories are made to facilitate a clinical interpretation of the results.