Juan David Ospina
email: <jdospina@gmail.com

Frédéric Commandeur

Richard Ríos

Gaël Dréan

Dr Juan Carlos Correa

Antoine Simon

Pascal Haigron

Renaud De Crevoisier

Oscar Acosta

Juan David Ospina

Crlcc Eug Ne Marquis

A tensor-based population value decomposition to explain

Introduction

Radiotherapy is one of the prescribed treatments for prostate cancer. Its objective is to deliver high doses of radiation to the tumour cells whereas sparing the neighbouring organs, often called organs at risk (OAR). The prediction of normal tissue complication has traditionally being addressed using dose volume histograms (DVH) 1 or dose surface histograms (DSH) which are reductions of the 3D [ ] dose distribution received by the organs. Although many studies have shown a correlation between dose, volume and rectal toxicity 2,3 , [ ] they lack spatial accuracy and are not able to correlate the treatment outcome with a specific dose pattern. One of the reasons is that normal tissue complication probability (NTCP) models mainly use the dosimetry information through the DVH, which does not preserve the spatial information. In radiotherapy is a well-known fact that similar DVHs may come from different dose distributions.

In the context of intensity modulate radiotherapy (IMRT) more adapted treatments could be proposed by including new constraints during the treatment planning step. Thus, the identification of anatomical regions correlated with toxicity is crucial to provide new recommendations for treatment planning. Some attempts have been made to introduce the notion of spatiality 4 , but still within the parametrizing the 3D patterns of dose. In this way, by selecting a limited set of predictive features, their method outperforms the classical models based on DVH/DSH. These approaches still lie in the reduction of feature dimensionality by fitting analytic functions to each dose shape. In 8 a principal components analysis (PCA) is developed. The PCA is applied on the registered dose matrices (expressed as row

[ ]
vectors) to obtain the principal modes of variation. Then, the modes of variation are used to classify patients according to a specific toxicity. The main issue of this work is the difficulty to interpret the results. Even if some principal components are related with a toxicity event, a spatial dose pattern is still hard to identify.

It would be desirable to have a method that allows for the identification of dose patterns. To accomplish this, we propose to generalise and use a population value decomposition (PVD) technique proposed in 9 for 2D image analysis. We propose a tensor-based

[ ]
generalisation of this PVD combined with a method to carry out statistical analysis that mimics the voxel-base morphometry (VBM) methods 10,11 . We show that the direct application of these methods, as proposed in 12 , is difficult mainly because of the normality

[ ] [ ]
assumptions for the test statistics. In the end we show that it is possible to exhibit dose patterns associated with some type of toxicity following prostate cancer radiotherapy treatment.

Tensor-based population value decomposition

The tensor-based population value decomposition presented herein is a generalization of the PVD presented in 9 . In this former

[ ]
reference, a PVD method is developed using a possibly-truncated singular value decomposition (SVD) along with a PCA over the individuals left-and right-singular eigen-vector matrices. This method enables compressing, thus leading to a more efficient ' representation of the subjects

Matrix-based population value decomposition

Consider a population consisting of 2D images, , 1 2 . The objective is to express as in eq. 1,

N Y ∈ ℝ i = , , , N … Y (1)
where and are common to the whole population and is specific to the -th subject. The steps depicted in 9 to obtain the

U W R i [ ]
individuals and the population matrices are as follows:

'

Obtain the SVD of each image .

Create the concatenated resulting matrices and . The rationale behind this procedure is that the space spanned by the columns of (or ) contains subject-specific left (or right) S D singular-vectors that, although they are not equal amongst them, they should be similar. Thus, applying PCA should allow for the calculation of the principal modes of variation of the left (or right) singular-vectors.

S = [S 1 | | … S ] D = [D 1 | | … D ]
In addition, a more efficient representation could be obtained as a previous step to the population analysis. This decomposition allows for two different compression levels:

At the individual s level some columns of and could be discarded regarding at some reconstruction error, leading to the As one can notice, the resulting individual specific matrices are not diagonal. This matrices can be regarded as projections onto a R subspace of , with and .

ℝ A ≤ F B ≤ G

Extension to 3D images

In the following approach the fact that a 3D image can be regarded as a third-order tensor is used. The SVD step of the 2D-PVD described before is emulated using a higher order singular value decomposition (HOSVD) 13 . The population consists of 3D images,

[ ] Y ∈ , 1, 2, , regarded as tensors. ℝ i = , N …
The procedure to obtain a 3D-PVD was implemented as follows:

Obtain the HOSVD of each image as in eq. ( 3),

where is the core tensor and are the corresponding factor matrices. For a 3D image

where .

As it can be seen from eq. ( 4), the image corresponding to the -th individual can be expressed in terms of a specific-individual core i tensor, namely , and population factor matrices. In addition, as in the 2D case, two compression levels are allowed:

Λ M
At the individual s level each image can be expressed as in eq. ( 5), '

where is a matrix consisting of the the first columns of , and is the truncated core tensor.

At the population level the matrices can be truncated by taking only the first columns leading to the reduced matrix . The

P Q P Q k
same reduction is possible when the PCA is carried out on .

Using both levels of compression the -th invididual s image can be approximated as in eq (6), i '

where is the -th individual s core tensor and , are the population factor matrices. i

' P Q 1 , … P Q M

Population analysis

Once each individual has been rewritten as in eq. ( 6) the statistical analyses can be carried out using the reduced individuals' representation , 1 . If the population can be split into two mutually exclusive groups, say and , an element-by-element

Λ i = , , N … G 1 G 2
comparison of the elements of could help to highlight the differences amongst both groups. This comparison can be done using a Λ two-sample t-test. The resulting p-values could be used then to decide which components are different amongst the two groups. Let Λ ( ) G 1 and be the typical core tensors of groups and , respectively, with components and . In addition, let be the

Λ ( ) G 2 G 1 G 2 Pv
resulting p-value tensor with entries . In this work the components of and are calculated as , 1, 2 if < pv

Λ ( ) G 1 Λ ( ) G 2 j = pv p
and otherwise, where is some specified threshold. It means that for the components where no statistical difference is p found the reconstruction is done using the information of the whole population whereas for those components where a statical difference is accepted the component reconstruction is done with the information coming from either one group or the other.

Using this approach the typical images for each group can be reconstructed as shown in eq. ( 7), ( 7) images consisted of 135 slices whose size in the axial plane was 512 512 pixels, with 1mm image resolution and 2mm slice thickness. The × median follow-up period was 38 months, with a minimum of 24 months for all patients. Rectal toxicity events were prospectively collected and scored according to the common terminology criteria for adverse events (CTCAE) version 3.0. The endpoint of the study was 2-year Grade 1 rectal bleeding, excluding acute toxicity. Patients with a history of hemorrhoids were not allowed to be scored as Grade 1 ≥ bleeding. In total, 12 patients presented at least a Grade 1 rectal bleeding event, which occurred between 6 and 24 months following the treatment.

Rectal bleeding study

Image processing

Patient s planned CT and dose distributions were non-rigidly registered with the demons algorithm 15 , on a single coordinate system was applied to the dose distributions in order to obtain the mapped dose distributions to a single template that were compared in this work. [START_REF] Kupchak | Experience-driven dose-volume histogram maps of NTCP risk as an aid for radiation treatment plan selection and optimization[END_REF] 5

The dose images were cropped in the axial plane to be of size 87 87 and only 51 slices were considered. Each dose was decomposed × using the HOSVD as shown in eq. ( 3) and no compression was allowed for at the individual or population levels as the interest was the precise identification of dose patterns.

Once all the patients core tensors were computed, an element-wise comparison was performed in the core tensor space using a ' two-sample t-test, considering the patients who presented at least a Grade 1 rectal bleeding event against to those who did not. A p-value lower than 0.05 was considered statistically significant. A typical core tensor for each group was calculated. Each element of the typical " " core tensors was computed as an average of the corresponding element across the population. For those elements that were found to be statistically different the average was computed using only the values of the corresponding group. For the remaining elements the average was computed using both groups values. The typical dose image for each group was reconstructed using the corresponding typical core " " " "

tensor and the population factor matrices. The difference amongst the two typical dose distributions was used to find regions correlated " "

with rectal bleeding Grade 1.

≥

Results

The Fig. 1(a) shows the difference between the typical bleeder individual and the typical non-bleeder individual obtained with the " " " "

proposed method. This result suggest an over irradiation (around 6 Gy) on patients presenting rectal bleeding Grade 1 in the zone of the ≥ posterior wall. For the sake of comparison in Fig. 1(b), the difference of the means of the bleeders doses and the non-bleeders doses is ' '

presented. The normality assumption was checked using a Shapiro-Wilk test in the core tensor space with a 95 confidence level. Only on % the 22 of the elements the normal assumption was rejected. This test was also performed in the native image space and the normality % assumption was rejected on the 67 of the voxels. Conversely, in the tensor case, no difference was detected on the elements where the % normality assumption was rejected. To illustrate the problem of normality assumptions, and then the validity of performing a two-sample t-test in the native space, Fig. 2 shows the histograms of the test statistics for both cases: the tensor case and the native image space. It is expected, if the normality assumption is valid, that the normalized histograms look like a standard normal distribution. This is the case for the tensor-based test statistics but not for those computed in the native image space.

Conclusions

This paper presented a method to perform population analysis for high order data ( 2) and its application to the study of rectal M > bleeding after prostate cancer radiotherapy. The presented method does not depend on the type of toxicity and it could be applied to find out other dose patterns correlated with different side effects, even in a different organ (such as the bladder). The main contribution of the tensor-based PVD was to enable the use of a classical statistic test to perform group comparison through the construction of the typical " "

individual for both groups. The obtained results should be carefully interpreted as the over irradiated zone depicted in Fig. 1(a) could depend on the selected template for registering all the images but also on the registration algorithm. Future work will include the segmentation of this zone and the study of the associated DVHs.

The direct application of the t-test over the dose images is not advised as normality assumption is rejected for most of the voxels. This is not the case in the core tensor space, however, highlighting the zones where both group are statistically different in the native image space is a subject of current research.

A future works will include as well a tensor-based PVD of each group and attempting to reveal the inter-group differences by the comparison of the population matrices. The effect of applying some level of compression at the individual s or population level and the ' validity of the method are also matter of future studies.

  Perform a PCA over to obtain the eigen-vector matrix of the covariance matrix such

  the matrix resulting of taking only the first columns of , is the analogous with respect to S L i L S D J i and consists only of the first columns and the first rows of . D V L i J i L J V At the population level it is possible to take only the first columns of (or the first columns of ), considering for example the A U B W explained percentage of variance, to build up the projection matrix (or ). The same reduction is possible when the PCA is performed U S L 1 … S L N ] D = [D J 1 … D J N ]

……

  Perform a PCA over each to obtain the projection matrices such as ( ) , 1, .S P S =P P S k = , MExpress each individual s image as in eq. (4),

DataA

  total of 63 patients treated for localized prostate cancer with IMRT were included in the study. The used treatment planning system was Pinnacle V7.4 (Philips Medical System, Madison, WI). The total prescribed dose was 46 Gy to the seminal vesicles delivered in 4.6 weeks, and 80 Gy to the prostate delivered in 8 weeks, with a standard fractionation of 2 Gy/fraction. The patient positioning, CT acquisition, volume delineations and dose constraints complied with the GETUG 06 recommendations as described in 14 . For the rectal [ ]wall, the constraints were maximal dose 76 Gy and V72 Gy 25 (the volume received by the 72 of the rectum wall volume). The CT

  by combining the CTs and organs delineations as explained in 8 . The patients CTs were first registered and then the same transformation [ ] '

 5 510 . Friston K , Holmes A , Worsley K , Poline J , Frith C , Frackowiak R . Statistical parametric maps in functional imaging: a general linear approach . Human brain mapping

Figure 2

Histograms of the test statistics used in the t-test.


Traian F Marinca

Ionel Chicinas

Olivier Isnard

Viorel Pop

Florin Popa

I Chicinaş

Synthesis, structural and magnetic characterization of nanocrystalline nickel ferrite-NiFe2O4 obtained by reactive milling

Keywords: nanocrystalline materials, magnetic materials, reactive milling, nickel ferrite

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

Introduction

Generally, all the soft magnetic ferrites are very important materials due to their magnetic and electric properties. The spinel ferrite suitable for many applications, such as: high frequency magnetic materials, recording media, ferrofluid technology, sensor technology, in microwave applications etc [1][2][3][4][5]. The soft magnetic ferrites are materials which crystallize in cubic spinel structure, spatial group Fd3m. The chemical formula of the spinel ferrite can be written as MeFe 2 O 4 , where Me is a metal or a group of metallic elements with two total valence. The spinel structure exhibits two kinds of crystal sites where the cations Me 2+ and Fe 3+ can be distributed: tetrahedral sites and octahedral sites. Typically two extreme type of spinel structure can be found : normal spinel and inverse spinel, but in practise mixed spinel structure are also observed. In normal spinel structure the tetrahedral sites are occupied by the Me 2+ and the octahedral sites are filled by the Fe 3+ cations. In inverse spinel structure the Fe 3+ cations are distributed between the tetrahedral and octahedral sites and the Me 2+ cations are found in octahedral sites. The mixed spinel structure has the cations distributed in a combined mode between the normal and inverse structure. The magnetic and electric properties of the soft ferrite are closely related to their intra-crystalline ordering of the cations [START_REF] Cullity | Introduction to Magnetic Materials[END_REF][START_REF] Goldman | Modern Ferrite Technology[END_REF].

Nickel ferrite (NiFe 2 O 4 ), like all soft magnetic ferrites, is the subject of intense researches in recent years. Nickel ferrite obtained by equilibrium methods has an inverse spinel structure with the Fe 3+ cations equally distributed between the two kinds of sites, Ni 2+ cations being located in octahedral sites and the magnetic ordering in the spinel structure of NiFe 2 O 4 is ferromagnetic. The iron cations from tetrahedral sites have the magnetic moments coupled antiparallel with the magnetic moments of nickel and iron from octahedral sites [START_REF] Cullity | Introduction to Magnetic Materials[END_REF]. It was highlighted that in the nanocrystalline/nanostructured state nickel ferrite has magnetic properties which differ from the polycrystalline state [START_REF] Kodama | [END_REF][9][10]. In order to produce nickel ferrite in nanocrystalline/nanostructured state several methods such as solgel [11,12], solvothermal [13], co-precipitation [14], bottom up [15], microemulsion [16] or hydrothermal [17] are used. Beside those methods the mechano-chemical synthesis route is a useful method for producing the nickel ferrite in nanocrystalline state. For producing the nanocrystalline nickel ferrite by mechano-chemical synthesis two ways are used: (i) in the first one, the ferrite is formed in the nanocrystalline state by solid state reaction between precursors during milling [9,[18][19][20]; (ii) in the second one, the polycrystalline nickel ferrite obtained by other methods is milled in order to reduce the grain size and to refine the structure [10,21,22]. Several studies where reported regarding the producing of nanocrystalline nickel ferrite by mechano-chemical route, but the results are very different depending on the experimental conditions. For example, the nickel ferrite phase was obtained after 8 hours [20] or else after 35 hours of milling [21]. For the sample milled 8 h a heat treatment was applied in order to increase the formation of nanoferrite [20]. Another study reported the obtaining of nickel ferrite in the nanocrystalline state by mechanical milling (with a crystallite size of 9 nm) from a mixture of oxides [2]. In the same article, it was observed that the magnetization is reduced by milling. A non-equilibrium distribution of cations and a spin canted effect have been reported by the same authors [2]. The cations inversion between the tetrahedral and octahedral sites into the spinel structure alongside of the spin canted effect and superparamagnetic behaviour was also reported by other authors [23]. The magnetization of nickel ferrite obtained in nanocrystalline state by mechanical milling is lower than the magnetization of nickel ferrite obtained by ceramic method [23]. Finally, it is worth to note that mechanical/reactive milling can be used to obtained nanostructures that have a nonequilibrium distribution of cations within the crystal structure.

The present paper presents the results regarding the synthesis of nanocrystalline nickel ferrite by reactive milling starting from an equimolar mixture of commercial oxides and the characterization of some structural and magnetic properties.

Experimental

A mixture of high purity commercial oxides powders (Alpha Aesar), nickel oxide (NiO -bunsenite) and iron oxide (α-Fe 2 O 3 -hematite) in molar ratio of 1:1 was used as starting sample. Stainless steel vials and balls were used for the milling process in high energy planetary mill Pulverisette 4 (Fritsch). The ball to powder mass ratio was 15: 1. The vial rotational speed (ω) and the disc rotational speed (Ω) were established at 800 rpm, respectively at -400 rpm. After milling, the samples were annealed in vacuum at a temperature of 350 °C for 4 hours in order to remove the internal stresses induced by milling and to improve the solid state reaction.

The formation of nickel ferrite by milling was investigated by X-ray diffraction (XRD) in the angular range 2θ = 25 -100° using a Siemens D5000 diffractometer which operates in reflection with CoK α radiation (λ = 1.7903 Å). The mean crystallite size was calculated using Scherrer's formula from full-width-at-half maximum (FWHM) of the diffraction peaks for the milled and annealed samples [24]. The lattice parameter was determined using Celref 3 program [START_REF] Laugier | CELREF V3. Developed at the Laboratoire des Matériaux et du Génie Physique[END_REF].

The morphology and chemical homogeneity of the particles were studied by scanning electron microscopy (SEM) using a JEOL -JSM 5600-LV microscope, equipped with energy dispersive spectrometer (Oxford Instruments, INCA 200 software). The particle size distribution has been determined using a Laser Particle Size Analyzer (Fritsch Analysette 22 -Nanotec) with an analysis field of 10 nm to 2000 μm.

Differential scanning calorimetry (DSC) investigations were performed using a alumina reference sample in NETZSCH-DSC 404S apparatus. The temperature range was between 25 and 1000 °C and a heating/cooling speed of 10 °C/min was used. The measurements were performed in Ar atmosphere in order to avoid the supplementary air oxygen.

The magnetization curves, M(H), were recorded at room temperature using the extraction sample method in a continuous magnetic field of up to 8 T [START_REF] Barlet | [END_REF]. The spontaneous magnetizations values have been derived from an extrapolation to zero field of the magnetization obtained in the magnetic field higher than 5.5 T.

Results and discussion

The X-ray diffraction patterns of the as-milled samples and those milled and subsequently annealed are shown in Fig. 1 (RMas-milled samples and RM + TT 350 ºC/4 h milled and subsequently annealed samples). For comparison, in the same figure is also shown the diffraction pattern for the starting mixture (SS). The Bragg reflections corresponding to the nickel ferrite (NiFe 2 O 4 ) are noticed in the diffraction pattern after 4 hours of milling, the NiFe 2 O 4 (220) and (511) peaks are clearly visible in the diffraction patterns at this time of milling. Those peaks are not the most intense peaks corresponding to the spinel nickel ferrite, the most intense peaks ((311) and ( 440)) are overlapped with the nickel oxides Bragg peaks. A broadening of the Bragg peaks of the precursors oxides, due to the decreasing of the crystallites sizes and to the second-order internal stresses induced by milling can be observed. A significant quantity of nickel ferrite was formed after 4 hours of milling. By increasing milling time, the Bragg peaks of the nickel ferrite are more clearly visible and those of the precursors oxides decrease. The most intense Bragg peak of the hematite phase, (200) peak, is almost invisible after 12 hours of milling. By increasing the milling time at 16 hours, the Bragg peaks corresponding to the nickel and iron oxides are no more present in the diffraction patterns; this indicates that the reaction between the precursors nickel and iron oxides is completed. The Bragg peaks intensity of the nickel ferrite phase does not follow the expectations reference JCPDS file . This highlights that nickel ferrite formed by milling has many inhomogeneities and is very disordered comparison to the ordered reference. The nickel ferrite structure is far from equilibrium, the spinel is rather a mixed one instead of a complete inverse. These results are in good agreement with earlier reported results [9]. By increasing the milling time, the system is getting far from equilibrium conditions and the nickel ferrite is forming continuously up to 30 hours of milling.

After only 4 hours of milling the (110) iron peak can be observed in the diffraction patterns. The iron peak intensity increases by increasing the milling time. This is a consequence of the powder contamination during milling with iron from the milling vial and balls [27]. After annealing, the Bragg peaks definition is better, the nickel ferrite peaks intensity are closer to the reference JCPDS file and the intensity of the peaks corresponding to the precursors oxides are reduced. This is due to the fact that the annealing removes the internal stresses induced by milling and increase the crystalline state of the synthesized NiFe 2 O 4 powder. For the milling time up to 12 hours, where the NiFe 2 O 4 is not completely formed, the annealing improves also the reaction between the milled and un-reacted precursors oxides [28]. The ratio of the most intense Bragg peaks intensities of the hematite and of the nickel ferrite for the 8 hours of milling is 0.59 and after annealing is decreasing at 0.32. So, the amount of nickel ferrite formed after annealing is significant. The evolution of the mean crystallite size and of the lattice parameter for the as-milled samples and respectively for milled and annealed samples versus milling time are shown in fig. 2. The crystallite size is typically in nanometer scale: 39 ± 2 nm for 12 hours of milling and decreases rapidly by increasing the milling time, up to 9 ± 2 nm after 30 hours of milling. The lattice parameters for the as-milled and for the milled and subsequently annealed samples increase by the increasing of the milling time. The lattice parameter of the as-milled samples is larger than for the milled and annealed samples, due to the internal stresses induced by milling and a non-equilibrium cations distribution in spinel structure. After annealing, the internal stresses are removed, a lattice relaxation occurs and it is possible to have a partial re-ordering of the metallic ions between tetrahedral and octahedral sites [29]. Whatever the milling time, the internal stresses and the ionic disorder lead to a larger lattice parameter of the milled ferrites than that of the nickel ferrite obtained by the ceramic method (a = 8.339 Å, JCPDS 10-0325 file). Another cause of the increased lattice parameter could be the higher iron content in the ferrite due to the iron contamination during milling. Indeed iron atoms could be incorporated into the nickel ferrite spinel structure. So, it is possible to have a nickel ferrite with more iron cations. Further investigations in order to clarify this aspect are in progress.

DSC curves recorded both on heating and cooling for the sample milled for 16 hours is shown in fig. 3. On the heating curve is remarked a large exothermic peak between 200 and 1000 °C. Released heat is associated with the lattice relaxation and with the cations reordering between tetrahedral and octahedral sites. Beside the large exothermic peak, on the heating curve an endothermic peak can be also remarked. This endothermic peak is associated with the Curie temperature of the nickel ferrite. The certainty that this peak is corresponding to the magnetic transition of nickel ferrite is given by the appearance of an exothermic peak at the same value (560 °C) on the cooling curve. The peak corresponding to the magnetic transition on the heating curve is larger, due to the distorted structure of the ferrite. Also the value of the magnetic transition temperature is lower compared with the values determined for the nickel ferrite synthesized by classical ceramic route (585 °C) [START_REF] Cullity | Introduction to Magnetic Materials[END_REF]. The structure of the obtained spinel ferrite and probably the impurities leads to a lower magnetic transition temperature. Typically the magnetic ordering in spinel nickel ferrite is ferrimagnetic [START_REF] Kodama | [END_REF]. The magnetic moments of the cations from tetrahedral sites are coupled antiparallel with the moments of the cations from octahedral sites. The magnetization is influenced by the disordered spinel structure produced by milling [9]. It may be noted that even for a magnetic field of 8 T the magnetization is not completely saturated (fig. 4). This behaviour could be attributed to the following phenomena: (i) the spin canted effect on the surface of the particles; (ii) a superparamagnetic behaviour due to very fine particles of ferrite obtained by milling and (iii) a combination of the two previously [START_REF] Kodama | [END_REF][9][10]. The M(H) curves of the milled samples are not saturated even after the annealing. Probably the spin canted effect is very strong and is not removed by the annealing or the un-saturated magnetization is caused in major part by the existence of very fine particles. The evolution of the spontaneous magnetization calculated for the starting sample (SS), as-milled samples and for the milled and annealed samples versus milling time is given in fig. 5. The formation of an important quantity of nickel ferrite after 4 hours of milling is confirmed by the magnetic analysis, the spontaneous magnetisation is 1.95 μ B /f.u. at this time of milling. For samples milled up to 12 hours, the spontaneous magnetization increases continuously by increasing milling time and it is reaching 10.84 μ B /f.u., as a consequence of the increased amount of ferrite phase in the material. After 16 hours of milling, when very probable the reaction between the precursors oxides is finished, the spontaneous magnetization decreases by increasing milling time. The spontaneous magnetization reaches 6.22 μ B /f.u. at the end of the milling process, after 30 hours. The magnetization values are lower in comparison with the value obtained for the nickel ferrite synthesized by ceramic route and in comparison with the calculated one [9,[START_REF] Cullity | Introduction to Magnetic Materials[END_REF]. These lower values could be attributed to the non-equilibrium distribution of the cations in spinel structure, structural disorder and spin canted effect [2,[START_REF] Kodama | [END_REF][9][10].

After annealing, the spontaneous magnetization increases in comparison with the magnetization of the as-milled samples. For the milling times where the NiFe 2 O 4 is not completely formed the magnetization increases due the increasing of the amount of spinel phase in material. For example, the spontaneous magnetization for the 4 hours milled sample increases at 3.42 μ B /f.u. after annealing. All the factors favouring the reduction of the magnetization (non-equilibrium cations distribution, structural defects) for the samples where the NiFe 2 O 4 was completely formed are removed or reduced by annealing, so the nickel ferrite is little closer to equilibrium and the crystallinity is increased [30]. The spontaneous magnetization is increasing after annealing with about 11 %, from 9.65 μ B /f.u. to 10.68 μ B /f.u. for the 16 hours milled sample. a b Fig. 6. SEM images of nickel ferrite powders obtained for the 16 h milled sample: a -500x; b -10,000x. The powder is formed by agglomeration of nanometric particles alongside by micrometric size particles.

The SEM studies show that the particles have a polyhedral shape, the powder consisting of micrometric particles and very fine nanometric agglomerated particles. In the Fig. 6 is presented a SEM images obtained for the 16 hours milled sample at a magnification of 500x and 10,000x. Very fine particles can be the single-domain particles with higher magnetization [START_REF] Cullity | Introduction to Magnetic Materials[END_REF]. The particle size distribution analyses confirmed the details obtained by SEM. There are particles of a few micrometers and particles whose size is in nanometric range. The median particle size D 50 determined for the 16 hours milled sample is 1. The EDX microanalyses confirmed the results obtained by XRD, which indicate a contamination of the samples with elemental iron. Whereas the starting sample EDX analysis indicated an iron/nickel ratio of 2:1, this ratio increases upon increasing the milling time. The amount of iron impurities is increasing upon increasing the milling time. For the 8 hours milled sample the Fe:Ni atomic ratio is increasing at 2.2-2.3:1 and after 16 hours of milling the atomic ratio of 2.45-2.6:1 is reached. This ratio becomes 2.6-2.8:1 for the sample milled for 30 hours. These ratios correspond to a powder contamination with iron between 3 and 10 percent. Small amount of elemental chromium was also detected by X-ray microanalyses. The chromium impurities are provided by the vial material which is consisting of stainless steel with 20 % of chromium. The contamination with chromium varying with the analysed micro-areas and increases by increasing the milling time. The amount of chromium detect in the sample is about 0.3 -0.4 % at for the sample milled 8 hours and reaches at one percent for the 30 hours milled samples.

The samples are chemically homogeneous and iron is uniformly distributed in the milled powder. In the fig. 8 is presented the EDX micro-analyses obtained for a micro-area from 16 hours milled sample. The chemical homogeneity is maintained upon increasing the milling time.

Conclusions

The nanocrystalline nickel ferrite was successfully obtained by reactive milling from a stoichiometric mixture of commercial nickel and iron oxides. After 16 hours of milling the reaction between precursor oxides is finished. The nickel ferrite is formed and support structural changes continuously up to 30 hours of milling. A very disordered nickel ferrite, with many defects was obtained by milling. The annealing improves the solid state reaction of the nickel ferrite formation and removes the internal stresses and the structural defects. The mean crystallite size decreases with increasing milling time, up to 9 ± 2 nm after 30 hours of milling. The lattice parameter decreases after annealing due to the spinel structure relaxation. DSC measurements revealed a large exothermic process associated with cations reordering in the crystalline structure. The magnetization is influenced by the disordered ferrite structure and by the internal stresses induced during milling. The formation of the spinel ferrite increases the magnetization, after the complete formation the magnetization decreases. The spontaneous magnetization increases after annealing for all the milling time, due to the reorganization of the cations into the spinel structure. A possible superparamagnetic component can be induced in the magnetization by very fine powder particles. The SEM images and particles size analysis confirmed the existence of very fine nickel ferrite particles. The processed powder was contaminated with iron and chromium from the vials and balls.

The reactive milling and subsequently annealing is an efficient route to synthesize nanocrystalline nickel ferrite powder. A little higher annealing temperature could lead to a higher magnetization without changing the nanocrystalline state of the material. Further experiments about this are in progress.

Fig. 1 .

 1 Fig. 1. X-ray diffraction patterns of the starting sample (ss), as-milled samples (4, 8, 12, 16, 20, 24 and 30 h of milling) and for the milled and subsequently annealed at 350 •C for 4 h samples. The position of the Bragg peak for each phase is indicated in the figure. For clarity, the diffraction patterns where vertically shifted.

Fig. 2 .

 2 Fig. 2. Evolution of the mean crystallite size and of the lattice parameter of the NiFe 2 O 4 for as-milled samples and respectively for milled and annealed samples versus milling time.

Fig 3 .

 3 Fig 3. DSC heating and cooling curves for the sample milled 16 hours. The peaks corresponding to the magnetic transition of the NiFe 2 O 4 on the heating and cooling curves are marked.

Fig. 4 .

 4 Fig. 4. Room temperature magnetization curves of the starting sample, milled for 4, 12, 16 and 30 hours samples and milled and subsequently annealed (350 °C/4 hours) samples.

Fig. 5 .

 5 Fig. 5. Spontaneous magnetisation of the starting sample, 4, 8, 12, 16, 20, 24 and 30 hours milled samples and of the milled and subsequently annealed (350 °C/4 hours) samples. For comparison some values from literature are given.

7 .

 7 1 microns and decreases by increasing the milling time up to 0.92 microns after 30 hours of milling. In the fig. 7a and fig. 7b are presented the particles size distributions for the 16 hours milled sample and 30 hours milled sample respectively. Particles size distribution of NiFe 2 O 4 powder milled obtained after 16 h of milling a) and respectively, 30 h of milling b).

Fig. 8

 8 Fig. 8 Distribution maps for Ni, Fe and O chemical elements obtained for the 16 h milled sample. The X-ray microanalysis was performed on the SEM image presented in the upper left side (10,000X).

Acknowledgments

This work was supported by CNCSIS -UEFISCSU, project number PNII -IDEI code 1519/2008.


Amal Mahboubi
email: amal.mahboubi@unicaen.frluc.brun

Luc Brun

Xavier Dupé

Object Classification Based On Graph Kernels

Keywords:

Automatic object recognition plays a central role in numerous applications, such as image retrieval and robot navigation. A now classical strategy consists to compute a bag of features within a sliding window and to compare this bag with precomputed models. One main drawback of this approach is the use of an unstructured bag of features which do not allow to take into account relationships which may be defined on structured objects. Graphs are natural data structures to model such relationships with nodes representing features and edges encoding relationships between them. However, usual distances between graphs such as the graph edit distance do not satisfy all the properties of a metric and classifiers defined on these distances are mainly restricted to the K nearest neighbors method. This article describes an image object classification method based on a definite positive graph kernel inducing a metric between graphs. This kernel may thus be combined with numerous classification algorithms.

INTRODUCTION

One now classical [START_REF] Harzallah | Combining efficient object localization and image classification[END_REF] strategy to recognize objects consists to compute a bag of features from a sliding window and to compare this bag to a model using scoring functions. Features of the bag may be deduced from boundary segments or SIFT points [START_REF] Harzallah | Combining efficient object localization and image classification[END_REF] defined at multiple scales.

However, such methods usually resume the image information to histograms hereby loosing any structural property of the object. Recent algorithms [START_REF] Xia | Pairwise Similarity Propagation based Graph Clustering for Scalable Object Indexing and Retrieval[END_REF] have proposed to overcome this drawback by encoding the objects by graphs using a heuristic graph matching algorithm to classify objects. Indeed, graphs provide a natural framework to capture the structure of objects. Classical distances between graphs are either based on the computation of the maximum common sub-graph [START_REF] Pelillo | Machine hierarchical structures using association graphs[END_REF] or on graph edit distances [START_REF] Riesen | Bipartie graph matching for computing the edit distance of graphs[END_REF] both problems being related and both nearly N P complete. The computation of a distance be-tween graphs using this approach requires thus either the use of an appropriate heuristic to compute efficiently an approximation of the distance on general graphs or a simplification of the graphs in order to apply more robust algorithms. Both solutions have been used. For example, Bunke [START_REF] Riesen | Bipartie graph matching for computing the edit distance of graphs[END_REF] proposed an efficient heuristic to approximate the edit distance between two graphs. On the other hand Pelillo [START_REF] Pelillo | Machine hierarchical structures using association graphs[END_REF] transforms shock graphs into shock tree in order to apply an efficient maximal clique detection algorithm.

Neglecting the complexity issue one additional drawback of graph matching or graph edit distances is the lack of properties of the resulting similarity or distance measures. Indeed, such measures do not usually provide a metric. As a consequence, usual graph distances cannot be readily combined with usual classification or regression algorithms coming from the statistical machine learning framework.

Definite positive graph kernels may be understood as a measure of similarity between graphs. Using Mercer's theorem such kernels are associated to a mapping function of graphs onto an Hilbert space such that the kernel between two graphs corresponds to the scalar product between their mapping within the Hilbert space. Using the famous "Kernel trick" many usual tools defined within the statistical machine learning framework may be used by using only kernels between graphs hereby avoiding the explicit computation of the mapping function which is usually unknown.

Several graph kernels have been proposed by Neuhauss

and Bunke [START_REF] Riesen | Graph Embedding in Vector Spaces by Means of Prototype Selection[END_REF] based on the edit distance. However such kernels are usually not definite positive. Dupé and Brun [START_REF] Dupé | Shape classification using a flexible graph kernel[END_REF] have proposed several kernels based on the notion of bag of paths. Outside the pattern recognition framework, many graph kernels have been proposed in the bioinformatic community [START_REF] Vert | Classification of Biological Sequences with Kernel Methods[END_REF]. Most of these kernels are based on kernels between bags of paths. This paper aims to define a kernel between graphs deduced from points of interest detected in two image to be compared. Within such a framework the main feature of both graphs is the point of interest together with its neighborhood. The notion of path is thus less relevant and we propose to associate to each graph a bag of vertices encoding the points of interest. The structure of the graph is captured by adding to each vertex the ordered sequence of its neighbors. Our graph construction scheme together with our new graph kernel are presented in Section . The proposed approach is validated through experiments in section .

GRAPH KERNEL CONSTRUCTION SCHEME 2.1. Graph Construction

As in the usual bag of word model, the first step of our method consists to extract a set of features from an image. The Scale Invariant Feature Transform (SIFT) [START_REF] Lowe | Distinctive Image Features from Scale-Invariant Keypoints[END_REF] provides a list of key points robust to changes in image scale, noise, illumination and local geometric distortion. Each key point is expressed by its x and y coordinates, the scale, the orientation and 128 numbers (in the range 0-255) representing the descriptor. Figure 1 illustrates a key point detection by the SIFT method on a face image. This unstructured set of points is then organized into a graph in order to take into account spatial relationships between key points. Given an image, we build a labeled graph G = (V, E, µ, F, w) where V, E and F respectively denote the set of vertices, edge and vertex's features while µ and w respectively correspond to a vertex labelling function from V to F and a weighing function from V to R+. Each key point is associated to a vertex v ∈ V , while its features µ(v) ∈ F are stored as a label of this vertex. The weighing function w encodes the relevance of each vertex within the graph. Within our framework, this weigh is defined for each vertex as the scale of its associated key point. This choice corresponds to the usual assumption that relevant and stable key points should be detected at large scales. Note however, that any relevant measure of the relevance of key points may be used as weigh. The set E of edges is defined from the k mutual nearest neighbors based on the x, y coordinates of key points: One edge (v, v ′ ) belongs to E if v belongs to the k nearest neighbors of v ′ while v ′ belongs to the k nearest neighbors of v. The degree of each vertex is thus bounded by k.

In order to take into account the local arrangement of vertices incident to a given vertex we explicitly encode for each vertex v the sequence of its neighbors encountered when turning counter-clockwise around v. The neighborhood N (v) of vertex v is thus defined as an ordered set (a sequence) (v 1 , . . . , v n ) of vertices. The first vertex of this sequence is arbitrary chosen as the upper right vertex (minimum of the sum of both coordinates).

Graph Kernel Construction

One classical way to define a kernel between graphs is to consider finite [START_REF] Dupé | Shape classification using a flexible graph kernel[END_REF] or infinite [START_REF] Kashima | Marginalized kernel Between Labeled Graphs[END_REF] bags of paths. Such kernels allow to consider global or semi global properties of graphs using walks traversing large part of them. However, using SIFT key point detector the information is mainly local and is thus essentially captured by graph's vertices and their immediate neighborhoods.

We thus define our kernel between graphs as a kernel between bags of oriented neighborhoods. The design of such a kernel implies the definition of a basic kernel between oriented neighborhoods. Considering the oriented neighborhoods N (u) and N (v) of two vertices u and v, with l u = |N (u)| and l v = |N (v)| , one straightforward solution consists to use a kernel between vectored data hereby considering each sequence as a vector:

K seq (u, v) = 0 if l u = l v lu i=1 K g (u i , v i ) otherwise (1)
where K g (u, v) is a RBF kernel between the features of input vertices:

K g (u, v) = e -d(µ(u),µ(v)) σ
, σ denoting a tuning parameter while d(., .) is the Euclidean distance, between the features µ(u) and µ(v) of vertices u and v. Sequences (u i ) i∈{1,...,lu} and (v i ) i∈{1,...,lv} denote respectively the oriented neighborhoods N (u) and N (v).

The kernel between oriented neighborhood defined by equation 1 compares the features of two oriented neighborhood of equal length. However due to acquisition noise or small changes between two images some SIFT points may be added or removed between two similar images. Such addition or removal of SIFT points may modify the neighbourhood of vertices belonging to the k nearest neighbors of the added or removed points. Using equation 1 such vertices become incomparable with their original version in the initial graph.

In order to overcome this drawback we introduce a rewriting rule on oriented neighborhood, which removes the vertex of lower weight of an oriented neighborhood. Given one vertex v, the rewriting of its oriented neighborhood N (v), denoted κ(v) is defined as:

κ(v) = (v 1 , ..., v i , ..., v lv ) (2)
where v i = argmin j∈{1,...,lv} w(v j ) is the neighbor of v with lowest weight. Note that, using a more rigorous notation, the rewriting of N (v) should be denoted κ(N (v)). We nonetheless use the simplified notation κ(v) in order to avoid cumbersome notations.

The function κ may be iterated thus leading to a sequence of successively simplified oriented neighborhoods: (κ 1 (v), κ 2 (v), ..., κ D (v)) where D is fixed to 50% of the maximal vertex's degree k (Section ). Each oriented neighborhood of such a sequence is obtained from the previous oriented neighborhood by removing one vertex. the dissimilarity with the initial neighborhood N (v) increases thus along the sequence. In order to quantify such a dissimilarity we associate a cumulative weigh to each element of this sequence of rewritings. The cumulative weigh of the i th rewriting of the initial oriented neighborhood N (v), denoted Cw(κ i (v)) corresponds to the sum of the weigh of the vertices removed to obtain κ i (v) from N (v) and is formally defined as:

Cw(v) = 0 Cw(κ i (v)) = w(v i ) + Cw(κ i-1 (v)) (3)
where v i is the vertex removed between κ i-1 (v) and κ i (v).

Based on these rewriting rules, we extend the initial vertex's kernel defined by equation 1 to a convolution kernel [START_REF] Haussler | Convolution Kernels on Discrete Structures[END_REF] between two sets of rewrote oriented neighborhoods:

K rewriting (u, v) = D i=1 D j=1 W (κ i (u), κ j (v)) * K seq (κ i (u), κ j (v)) (4)
where the kernel W penalizes costly rewritings corresponding to the removal of important key points. This kernel is defined as a kernel between the cumulative weighs of both rewritings:

W (κ i (u), κ j (v)) = e -Cw(κ i (u))+Cw(κ j (v)) σ ′ (5)
where σ ′ is a tuning variable. Our final kernel between vertices combines K rewriting with one kernel between the features of the two central vertices defining both oriented neighborhoods:

K(u, v) = K g (u, v)K rewriting (u, v) (6)
Our final kernel between both graphs corresponds to a kernel between both bags of oriented neighborhoods using equation 6. One usual design of such a kernel [START_REF] Dupé | Shape classification using a flexible graph kernel[END_REF] would consist to perform a sum of kernels between all couples of oriented neighborhoods belonging to both bags. However, such a solution do not take into account the relevance of vertices within both graphs. In order to take such a relevance into account we introduce a function ϕ from V to R+ which increases with respect to vertex's weigh and is formally defined as:

ϕ(v) = e - 1 σ ′ (1+w(v)) .
Our final graph kernel is thus defined as a weighed sum of kernels between oriented neighborhoods:

K graph (G 1 , G 2 ) = u∈V1 v∈V2 ϕ(u)ϕ(v)K(u, v) (7)

EXPERIMENTS

Dataset

We choose to test our kernel against a database of 41 faces composed of 11 females and 30 males, extracted from the FERET face database 1 . Each of the 41 subject of the database is associated with 5 images: one regular facial expression, an alternative facial expression, two half left and half right rotations of -67.5 degrees and one right profile with a rotation of +90 degrees. The image is fixed to 256 × 384 pixels. Some images from the database are displayed in Figure 2.

Regular expr.

Alternative expr. Profile image

Experimental Setup

Our first experiment illustrates the importance of the rewriting step within our kernel computation scheme.

To this end, we randomly select 5 images among our database and corrupt them by an additive Gaussian noise with increasing standard deviation σ. The similarity between each original image and its corrupted version is measured both with our kernel K graph (equation 7) and by a kernel without any rewritings denoted K graphSeq . This last kernel is defined from equation 7 by setting

K(u, v) = K g (u, v).K seq (u, v) (equation 1).
In order to compare both kernels we normalize them using the following formula [START_REF] Haussler | Convolution Kernels on Discrete Structures[END_REF]:

K(G, G ′ ) = K(G, G ′ ) K(G, G)K(G ′ , G ′ ) (8)
where K represent the normalized version of a kernel K. Using such a normalisation step, both kernels Kgraph and KgraphSeq vary between 0 and 1. the standard deviations of kernel values computed on the 5 images for each noise levels. As shown by Fig. 3, the curve representing KgraphSeq quickly drops to 0. This results means that even for low noise level, any SIFT point in a corrupted image has in its neighborhood at least one added or removed point compared to the same SIFT point in the original image. Using equation 1, the similarity between the oriented neighborhoods of both SIFT points is equal to 0. On the other hand, the curve associated to Kgraph quickly drops to 0.2 for a value of σ equal to 1, and then softly decreases as the noise level increases hereby showing the resistance to noise bring by our rewriting step.

Our second experiment investigates the relationships between the number of rewritings D and the maximal size of the neighborhood k. To this end, we randomly selected 4 images belonging to different classes within the FERET database. We then iteratively removed each of the 4 selected images from the FERET database and computed its 4 nearest neighbors among the 201 remaining images. For a given choice of k and D, if our kernel correctly captures the class of each image, the 4 closest neighbors of each query should correspond to the remaining images of its class. We thus computed for each of the 4 query images the number of images belonging to its class among its 4 closest neighbors within the database. This score computed on each query image is sum up in order to obtain a global result for each couple (k, D) which varies between 0 and 16. Table 1 summarizes this experiment for different choice of k and D. The maximal value of each line is highlight by bold fonts. One can note that for each value of k the score passes through a single maximum. This phenomenon may be interpreted as follows: For low value of D the number of rewritings do not compensate small image perturbations due to noise and view change. Conversely, for high,value of D, our rewriting step may remove important points with slightly similar weigh in both neighborhoods. Such removals may lead to different neighborhoods with a low similarity according to K seq (equation 1). The weighted sum of this low values of K seq leads to a low overall value of K rewriting (equation 4). The optimal value of 1 2 3 4 5 1 3 1 1 0 0 2 0 5 0 0 0 3 0 0 5 0 0 4 0 1 2 2 0 5 0 0 2 0 3 D is equal to 50% of k for k = 10 and k = 15. This optimal value remains equal to 5 for k = 7 and 13 hence representing respectively 71% and 38% of k. Therefore, as expected, the heuristic setting D = K 2 made in Section does not provide the optimal value of D in all cases. This choice provides however results close to the optimal solution and may be used as an alternative to an explicit estimation of the optimal value of D using a cross validation step.

Object Classification

Given a query image of an object of interest, object classification consists to retrieve its class from a set of objects already classified called a training set. Our graph kernel has been tested within this framework on 5 classes of the FERET database using the classical KNN classification algorithm. Figure 4 represents the confusion matrices obtained on these five classes for K = 1 and K = 4. As illustrated by Figure 4 both matrices have only few out diagonal values, the best results being obtained for K = 1. An example of the 5 nearest neighbors of a query image is represented in Fig. 5. Three of the four images belonging to the same class as the query and encoding the same person in different views are correctly retrieved. The missing image corresponds to the profile view which induces a very different graph.

CONCLUSION

In this study, a graph kernel method for object retrieval has been proposed with an application to face recognition. We defined a graph kernel based on bags of oriented

Figure 1 .

 1 Figure 1. SIFT Feature Points

Figure 2 .

 2 Figure 2. FERET Database

Fig. 3 Figure 3 .

 33 Figure 3. Normalized Kernels Kgraph and KgraphSeq .

Figure 4 .

 4 Figure 4. Confusion Matrices Using The KNN Classification Algorithm, 5 Classes Of The FERET Database And Two Values Of K.

Figure 5 .

 5 Figure 5. Ranked Five Closest Images From the Query

Table 1 . Classification Scores Obtained On The FERET Database For Various Values Of k

 1 And D.

	D=1 D=3 D=4 D=5 D=6 D=7 D=12
	k=7 3	5	5	7	6	-	-
	k=10 3	5	5	6	5	6	-
	k=13 3	5	5	7	5	6	6
	k=15 6	5	5	5	5	7	6


Jun Wang

Shijie Wang

Yang Chen

Jiasong Wu

Jean-Louis Coatrieux

Limin Luo

Keywords: computed tomography, metal artifact reduction, image fusion 5>!

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

Metal artifact reduction in CT using fusion based prior

image.

projections [START_REF] Bal | Metal artifact reduction in CT using tissue-class modeling and adaptive prefiltering[END_REF][START_REF] Müller | Spurious structures created by interpolation-based CT metal artifact reduction[END_REF][18][START_REF] Wu | Metal artifact reduction algorithm based on model images and spatial information[END_REF][START_REF] Karimi | Segmentation of artifacts and anatomy in CT metal artifact reduction[END_REF][START_REF] Boas | Evaluation of two iterative techniques for reducing metal artifacts in computed tomography[END_REF][START_REF] Jeong | Reduction of artifacts due to multiple metallic objects in computed tomography[END_REF][START_REF] Prell | A novel forward projection-based metal artifact reduction method for flat-detector 56?! computed tomography[END_REF][START_REF] Prell | Metal artifact reduction for clipping and coiling in interventional C-arm CT[END_REF] . The first step is to determine the metal trace occupied by the corrupted data in the sinogram space. To >?! achieve this, a re-projection method is often used. An image is first reconstructed from the original sinogram. Then, the metal parts are identified in the image domain by thresholding. Forward projection of the metal-only image yields the metal-only sinogram, in which the nonzero pixels define the metal trace. A correction can then be performed to estimate the surrogate data in the metal trace. Correction methods can be divided into two groups: interpolation-based methods [6][7][START_REF] Gu | A method based on interpolation for metal artifacts reduction in CT images[END_REF][START_REF] Gu | X-ray CT metal artifacts reduction through curvature based sinogram inpainting[END_REF][START_REF] Bertram | Directional view interpolation for compensation of sparse angular sampling in 5">! cone-beam CT[END_REF][START_REF] Yan | Projection correlation based view interpolation for cone beam CT: primary fluence restoration in scatter measurement with a moving beam stop array[END_REF][START_REF] Wouter | Development and validation of segmentation and interpolation techniques in sinogram for metal artifact suppression in CT[END_REF][13][START_REF] Raoul | Metal artifact reduction for CT: Development, implementation, and clinical comparison of a generic and a scanner-specific technique[END_REF][START_REF] Kratz | A fully 3D approach for metal artifact reduction in computed tomography[END_REF] and prior-based methods [START_REF] Bal | Metal artifact reduction in CT using tissue-class modeling and adaptive prefiltering[END_REF][START_REF] Müller | Spurious structures created by interpolation-based CT metal artifact reduction[END_REF][18][START_REF] Wu | Metal artifact reduction algorithm based on model images and spatial information[END_REF][START_REF] Karimi | Segmentation of artifacts and anatomy in CT metal artifact reduction[END_REF][START_REF] Boas | Evaluation of two iterative techniques for reducing metal artifacts in computed tomography[END_REF][START_REF] Jeong | Reduction of artifacts due to multiple metallic objects in computed tomography[END_REF][START_REF] Prell | A novel forward projection-based metal artifact reduction method for flat-detector 56?! computed tomography[END_REF][START_REF] Prell | Metal artifact reduction for clipping and coiling in interventional C-arm CT[END_REF] . >>! Interpolation-based methods estimate the surrogate data using interpolation type techniques. Linear interpolation (LI) is the simplest way to do that 6,[START_REF] Gu | A method based on interpolation for metal artifacts reduction in CT images[END_REF] . We refer them to LI-MAR in this work. In LI-MAR, the surrogate data are estimated via linear interpolation using the uncorrupted data on the two sides of the metal trace along the channel ray??? direction.

Interpolation-based methods can effectively reduce the metal artifacts. However, due to the unavoidable loss of structure information in the metal trace after interpolation, they often introduce new artifacts near the metallic objects. @?! To overcome the drawbacks of interpolation-based methods, prior-based methods estimate the surrogate data using forward projections from a prior image. The structure information from the prior image is brought back to the metal trace and this way, better estimates can be obtained. Most existing methods, based on k-means clustering [START_REF] Bal | Metal artifact reduction in CT using tissue-class modeling and adaptive prefiltering[END_REF][START_REF] Wu | Metal artifact reduction algorithm based on model images and spatial information[END_REF] or multithresholding [START_REF] Müller | Spurious structures created by interpolation-based CT metal artifact reduction[END_REF]18,[START_REF] Prell | A novel forward projection-based metal artifact reduction method for flat-detector 56?! computed tomography[END_REF][START_REF] Prell | Metal artifact reduction for clipping and coiling in interventional C-arm CT[END_REF] , build the prior images by segmenting the originally reconstructed images (the uncorrected images) or pre-corrected images (via interpolation-based MAR) into different tissues such as air, soft tissue, and bone. @>!

We will call them TP-MAR methods since all are based on thresholding the prior images. One main limitation of TP-MAR methods is their sensitivity to misclassification which can lead to some tissue structures lost. In Ref. 21, a metal deletion technique (MDT) was proposed that used a LI-MAR corrected image as an initial prior image. Then, the prior image was iteratively improved using the FBP reconstruction combined with an edge-preserving blur filter. Since MDT was still based on the interpolation-based MAR, tissue structure loss near the metal objects was often observed.

A?!

Another limitation of MDT is the high computational cost, about 19 times slower than FBP. Recently, a method based on computer vision techniques was developed [START_REF] Karimi | Segmentation of artifacts and anatomy in CT metal artifact reduction[END_REF] . This method first discriminated the artifacts from the real tissue ! structures in the uncorrected image using shape and intensity information of metal artifacts. Then, the prior image was generated by reassigning the CT values within the artifact region to the surrounding tissue values. The experiments on head CT scans reported in Ref. 20 has shown that this method could generate more accurate prior images than TP-MAR A>! methods, and could provide better image quality than both TP-MAR and interpolation-based MAR. Even though, its performance and robustness, in the presence of different metal implants or more pieces of metals, still needs a further verification, as pointed out by the authors.

667&!B0C:;

The quality of the prior image plays a crucial role in the performance of the prior-based MAR methods. To B?! overcome the misclassification resulting from segmentation in TP-MAR methods, we choose to build the prior image via the fusion of a pre-corrected image and the uncorrected FBP image. The main idea behind this new method (called fusion prior-based MAR or FP-MAR in the following) is to jointly benefit from the well preserved structures in the uncorrected FBP image and the suppressed artifacts in the pre-corrected one. We can see in Fig. 1

667=7&D'"E*-''"*#(-.&

The pre-corrected image is here obtained by applying LI-MAR and next, an edge-preserving blur filter. The uncorrected FBP image is first reconstructed from the original sinogram. Then, the re-projection method (see Sec. I. B) is used to identify the metal trace. The metal parts are segmented by applying a simple threshold, set to 2000 Hounsfield unit C>! (HU) . For dental filling cases, 3000 HU is used because human teethes usually have CT values above 2000 HU. After determining the metal trace, the linear interpolation technique proposed in Ref. 8 is used for sinogram completion. Fig. 1(a) shows the uncorrected FBP image of a simulated hip phantom with bilateral prostheses which suffers a prominent dark streak artifact between the two prostheses. Fig. 1(b) shows the result of LI-MAR, and it can be seen that the streak artifact is significantly reduced but simultaneously newly artifacts are introduced and serious loss of the bone structures "??! near the prostheses can be noticed.

An edge-preserving blur filter is then applied to suppress the newly introduced artifacts in the LI-MAR image. Each value ! !" !!! !! in the filtered image is the average value of pixels within a local square region around the pixel !!! !! in the LI-MAR image. To preserve tissue edges, only those pixels with CT numbers similar to the value of the center pixel are included in the average:

"?>! !!!!!!!!!!!!! !" !!! !! ! !!!!!!!!!! ! !!!! ! !!!! ! , ( 1
)
where ! is the number of "similar pixels" in the square region of size ! with the indices!!! and q running from ! ! to !. !!! ! !! ! ! !! is defined as:

!!!!!!!!!!!! ! ! !! ! ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! !" ! ! !! ! ! ! ! ! !" !! ! ! ! ! !" ! ! !! ! ! ! !!!!!!!! ! !" ! ! !! ! ! ! ! ! !" !! ! ! ! , ( 2
)
where ! !" denotes the CT values in the LI-MAR image. ! is a user-defined threshold for the identification of the ""?! similar pixels. ! and ! should be properly selected to achieve artifact suppression without blurring the main edges between air, soft tissue, and bones. Our experiments demonstrate that!!! ! !""!!" !"#!! ! !! can lead to good results for the all cases studied in this work. Fig. 1(c) illustrates the filtered result of the LI-MAR image (henceforth referred to as the pre-corrected image) where the newly introduced artifacts are well suppressed.

! @! ! 667@7&65$2"&),1(-.& "">!
This step aimed at building the prior image from the fusion of the pre-corrected image and the uncorrected FBP image with metal parts removed. Fig. 1(d) shows the resulting metal-removed FBP image, in which the metal parts are filled by the CT numbers taken from the metal regions in the pre-corrected image. This metal-removed FBP image is then subtracted from the pre-corrected image to generate the difference image ! depicted in Fig. 1(e). The positive values with large magnitudes in ! correspond to the dark streak between the two prostheses, whereas the negative "*?! values mainly represent the damaged bone structures. The prior image, ! !"#$" , is then built via the following pixel-wise fusion formula:

! !"#$" !!! !! ! ! !! ! !! !"# !!! !! ! ! ! ! !! ! !! !"#$ !!! !!, (3)
where ! !"# and ! !"#$ denote the metal-removed uncorrected FBP image and the pre-corrected image, respectively.

! !! ! is the pixel-wise weighting factor in the interval !! ! , ]0,1] ?? and is calculated through Eqs. ( 4) and ( 5):

"*>! !!!!!!!!! !! ! ! ! !! ! !"#$ !!! ! ! !, (4)
!!!!!!!!! !"#$ !!! !! ! ! ! !!! !!!! !"# ! !"# !!!! !"# , (5)
where

!! ! ! ! ! and ! ! !. ! !"#$ is
!! ! !!"!!"#$!!"# small values of !!! !"#$ !! ! !! in order to get ! !"#$" !!! !! close
to the intensity in the metal-removed uncorrected FBP image, thus preserving tissue structures. "5>! Parameters ! and ! should be suitably set to ensure both effective tissue preservation and artifact suppression for the prior image. Figs. 3(a)-3(d) show that a small ! may lead to severe residual dark artifacts, and the fusion is no longer sensitive to ! when ! is larger than 10 in this case (little difference is observed when increasing !!!from 10 to 40). Setting parameter ! is application-dependent. A small !! value would generate a fusion image similar to the pre-corrected image, while a large ! would make the result close to the metal-removed uncorrected FBP image (see Fig.

"6?! 3, bottom row). Thus, a small ! value (for example, 0.1) should be retained for cases where small and regular metal objects are present (such situations are often the ideal applications for LI-MAR), while it should be increased (for example, to 0.45) for irregular and large metal objects to overcome the possible serious loss of tissue structures near the metal objects. Simulations were performed to evaluate the proposed FP-MAR method. As shown in Fig. 4, we generated two phantoms with metal artifacts produced from bilateral hip titanium prostheses and from pedicle screws, respectively.

The original true phantoms without metallic objects are also displayed in Fig. 4 to serve as reference. With 2D fan beam scans assumed, the simulations modeled the CT scanning data using the CatSim model [START_REF] De Man | CATSIM: A new computer assisted tomography simulation environment[END_REF] which incorporates polychromaticity, realistic quantum and electronic noise models, non-linear partial volume effect, and scatter. The main "@?! simulation parameters were 70 kV, 10 8 initial photons per ray, 10 3 scatter photons per detector element, 768 channels, and 720 views per rotation. 10 rays were simulated per detector to preserve the finite width of beams.

As presented in Fig. 6 and Fig. 7, seven clinical cases were selected to test the capability of the proposed FP-MAR method in coping with different kinds of metal artifacts in standard clinic. All the original uncorrected images, "@>! (http://www.revisionrads.com/cgi-bin/MDT_cgi). To be more precise, the metal artifacts in the three sets of clinical images in Fig. 6 are respectively caused by a deep brain stimulator (Case 1), bilateral hip prostheses (Case 2), and small "A?! metal implants (Case 3). Compared to Fig. 6, Fig. 7 includes four more challenging cases with multiple metal implants of large size and irregular shapes. The metal artifacts of these four cases are respectively produced from unilateral hip prostheses (Case 4), sternum fixation screws (Case 5), pedicle screws (Case 6), and dental fillings (Case 7). For each case, the original sinogram was simulated by forward projecting the tested clinical images. The resolution of the CT image is !"#!!"# pixels and the simulated sinogram is of !"#!!"## pixels (3600 views over 360°).

"A>! ! C! !
In the implementation of the proposed FP-MAR method, ! ! !" and!!! ! !!!" were used in simulations. For the clinical image experiments, ! ! !" was used in all cases; t was set to 0.1 for the three cases reported in Fig. 6 and to 0.45 for the four challenging cases illustrated in Fig. 7. For comparison, the LI-MAR and TP-MAR methods were also implemented. We applied LI-MAR by using the pre-correction step procedure above (Sec. II. A)D! To generate the TP-MAR prior image, the uncorrected FBP image was first smoothed using the edge-preserving blur filter described in "B?! Sec. II.! ED Then, the filtered result was segmented into air, soft tissue, bone, and metals using the k-means clustering described in Ref. 16 with the initial clustering center -950 HU, 200 HU, 750 HU, and 5000 HU. After segmentation, the air regions were set to -1000 HU, the soft tissue and metal parts to 0 HU while bone pixels kept their original intensities as suggested in Ref. 18. Then, the data replacement for the metal trace was completed by the above difference-based interpolation method described in Sec. II.! 2. Since the interpolation technique is the same for all methods, the prior "B>! image determines the improvement. In addition, we collected the clinical MDT results from the REVISION RADIOLOGY group for comparison. All the CT images were processed under a PC workstation (Intel Core 2 Duo CPU 2.93 GHz processor and 2048 Mb RAM) under MATLAB R2009a version 7.8.0 (The MathWorks Inc., Natick, MA).

The fan-beam transform and the inverse fan-beam transform were used for the forward projection and filtered back projection, respectively. "C?! 6H7&9BA<F0A& Fig. 5 shows the correction results for the two simulations. We can see that the uncorrected images suffer from severe beam hardening artifacts between the metal objects and all the LI-MAR, TP-MAR, and FP-MAR methods lead to significant reduction of metal artifacts. However, in the LI-MAR results, some new artifacts are introduced and some bone structures (indicated by black arrows) are missing near the metallic objects. When compared to LI-MAR, "C>! TP-MAR yields better preserved bone structures, but some residual artifacts (indicated by white arrows) are also present due to segmentation errors in the prior image. We can see in Fig. 5 that the proposed FP-MAR method leads to the best image quality with effective artifact suppression and bone structure preservation. The most challenging cases are illustrated in Fig. 7. Here too, the LI-MAR method reduces metal artifacts at the cost of introducing new artifacts and losing tissue structures near the metal objects. The quality of LI-MAR images is even lower than the quality of the uncorrected ones. This stems from the fact that large metallic objects lead to large and (C=200 HU/W=2000 HU). Case 6: (n=10/t=0.45) and (C=200 HU/W=800 HU). Case 7: (n=10/t=0.45) and (C=400 HU/W=2000 HU). metal trace region, resulting in serious loss of structure information in the metal trace after linear interpolation. When compared to LI-MAR, TP-MAR works better in both preventing the introduction of new artifacts and preserving tissue structures near the metals. But there are still some dark or bright artifacts remaining around the metal objects in **>! TP-MAR processed images (see the regions indicated by dotted arrows) due to a misclassification in building the thresholded prior image. MDT leads to effective artifact suppression but tissue structures are missing near the metal objects (these regions are defined by dotted arrows). Obviously, the proposed FP-MAR performs better than LI-MAR and TP-MAR in suppressing artifacts and also to MDT in restoring those tissues close to metallic objects. H7&;6A/<AA6:8& *5?! ! "5! ! TP-MAR methods build the prior images by segmenting the uncorrected FBP image into air, soft tissue, and bone.

Streak artifacts can be as dark as air or as bright as bone, therefore segmentation errors cannot be avoided (refer to Fig. 8, second column) and they most often lead to residual artifacts in the final corrected images. Generally, less segmentation errors can be expected when segmenting an image that has been pre-corrected by means of interpolation-based MAR. Unfortunately, the quality of pre-corrected images may not be guaranteed, particularly when *5>! large metal objects are present. For example, the normalized metal artifact reduction 18 (NMAR) recommended the LI-MAR as the pre-correction method, but the image pre-corrected with LI-MAR can still suffer heavy residual artifacts and tissue structure missing as shown in Fig. 7. Using poor pre-corrected images undoubtedly yields poor prior images (see Fig. corrected image, the contribution of the two images being weighted according to their difference information. In Fig. 8, ! we can see that the proposed FP-MAR prior images, compared to TP-MAR priors, better preserve tissue information and significantly reduce residual artifacts. This is the main reason why FP-MAR outperforms TP-MAR. The proposed image fusion technique can also be used to provide a better pre-corrected image for existing TP-MAR methods. Figure 9 *>?! illustrates two thresholded prior images respectively segmented from the LI-MAR corrected image (middle) and from the proposed fusion image (right). The LI-MAR based prior image suffers from obvious residual artifacts and bone structure loss (see the dotted arrows). In contrast, the proposed fusion image leads to more accurate segmentation with almost no residual artifacts and much better preservation of bone structures.

The main limitation of FP-MAR comes from the fact that some bright artifacts around the metal objects can be *>>! wrongly incorporated into the prior image as bone structures (Fig. 8, black dotted arrows) and so, can lead to residual artifacts in the final corrected image. Decreasing the parameter ! can alleviate this problem by making the prior image closer to the filtered LI-MAR image. However, this operation may result in tissue lost if some original structures fail to be preserved in the filtered LI-MAR image (see the fusion image in Fig. 3 with ! ! !!!). Future work is therefore needed to address this limitation. The computer vision techniques, as mentioned in Ref. 20, may have the potential for *@?! solving this problem.

Several other limits of the FP-MAR method should also be addressed. First, the simple thresholding used to identify the metallic parts might lead to inaccurate identification of the metal trace, and this may introduce some fine streak algorithms must therefore be considered in the future to improve the segmentation. Second, since the fusion is performed on pre-corrected images using interpolation-based MAR method, improvements can be anticipated from more advanced interpolations [START_REF] Bertram | Directional view interpolation for compensation of sparse angular sampling in 5">! cone-beam CT[END_REF][START_REF] Yan | Projection correlation based view interpolation for cone beam CT: primary fluence restoration in scatter measurement with a moving beam stop array[END_REF][START_REF] Kratz | A fully 3D approach for metal artifact reduction in computed tomography[END_REF] . Third, if the used edge-preserving blur filter suppresses artifacts, normal tissues are *A?! simultaneously blurred (Fig. 1(c)), thus lowering the quality of the prior image. To further improve the quality of prior images, some advanced edge-preserving filters with artifact suppression property [START_REF] Bal | Metal artifact reduction in CT using tissue-class modeling and adaptive prefiltering[END_REF][START_REF] Chen | Thoracic low-dose CT image processing using an artifact-suppressed 5@?! large-scale nonlocal means[END_REF] should be considered. Fourth, several parameters need to be manually set in the current method (e,g. the threshold in metal segmentation, the!!! !"#!! in Eqs. (1) and (2), the ! and ! in Eq. ( 4)). Among all the parameters, ! requires a special attention in order to improve the prior image. Finally, due to the additional forward projection, FBP reconstruction, filtering, and fusion *A>! steps, extra computation is required for the proposed FP-MAR method: the whole execution time is about 216 seconds, which is approximately 4.5 times slower than one FBP reconstruction. Nevertheless, the overall algorithm can be accelerated by using GPU-based parallelization 34 to fulfill the clinical requirements (no more than 2 seconds per 2D image).

  FIG. 1. The intermediate outputs in generating a prior image of a simulated hip phantom with bilateral prostheses. The uncorrected FBP image, the LI-MAR image, and the pre-corrected image are shown in (a), (b), and (c), respectively. In (d), the uncorrected FBP image with metal parts removed. The difference image [subtract (d) from (c)], and the prior image obtained by fusion of (c) and (d), are depicted in (e)

  FIG.2. Weighting function corresponding to Eq. (4) with parameters n=10 and t=0.45.

  FIG. 3. Illustration of the fusion result for different settings of ! and !. Window settings: C=400 HU/W=1500 HU.

FIG. 4 .

 4 FIG. 4. The two simulated phantoms. The left and right sides depict respectively the phantoms with and without inserted metallic objects. Window settings: C=100 HU/W=700 HU.

!

Fig. 6

 6 Fig.6illustrates the correction results for Case 1, Case 2, and Case 3. Severe streak artifacts tangent to the metallic

  FIG. 7. Another four sets of correction results for Case 4-Case 7. Case 4: (n=10/t=0.45) and (C=100 HU/W=800 HU). Case 5: (n=10/t=0.45)

  FIG. 8. Prior image illustration of TP-MAR and FP-MAR methods for the experiments with the hip phantom (clinical Case 4 and Case 7). TP-MAR prior images suffer from obvious segmentation errors, as indicated by the white dotted arrows. The proposed FP-MAR prior

FIG. 9 .

 9 FIG. 9. Two thresholded prior images corresponding to Case 5 respectively segmented from the LI-MAR corrected image (middle) and from

H67&/: 8 /

 8 F<A6:8& *B?!The proposed FP-MAR approach described in this paper demonstrates its clinical potential by providing better corrected images than LI-MAR, TP-MAR and MDT methods in a wide range of cases with single or multiple pieces of metals in various shapes and sizes. Several improvements are under progress to face the limitations discussed above. In particular, we generate artificial sinogram from the tested clinical image as a base for correction. The fan-beam transform and inverse fan-beam transform have a low-pass filtering effect on the final corrected image, and this will *B>! lead to lowered spatial resolution. Future work is therefore needed to verify FP-MAR on real measured projections. Je ne comprends pas cette partie! Il faudrait la reformuler et ce qui serait bien qu'elle degage quelques idées d'amélioration. =/I8:JFB;G!B80A& This work was supported by the National Basic Research Program of China, No.2010CB732503. We are very *C?! "A! !