Lotfi Berguiga

Cyril Lupi

Nicolas Hyon

Laurent Salomon

Normand " Guy

Roger Le Ny

Eric Tanguy

Dominique Leduc

Monique Auvray

Jean Claude Bizeul'

Alain Gouronnec'

Frédérique De Fornel

Christian Boisrobert°

Michel Poulain

Marcel Poulain

Accelerated aging and life tests of optical fibers

Keywords: silica, aging, reliability, interferometry, microscopy

INTRODUCTION

The reliability of the network optical fibers is a critical issue for telecommunications.New investigations methods have been developed within a cooperative program supported by France Telecom. They include low coherence interferometry and optical near field microscopy in association with classical analysis tools such as Scanning Electron Microscopy and Shear Force Microscopy. One aim of the study is the localization and the characterization of the defects from which failure originates. The aging effect has been investigated in silica fibers immersed in desionized water at 65 oc and 85 "C for different times: from 3 to 12 months 1. EXPERIMENTAL SET-UP 1.1 Near field optical microscopy Light emitted by objects consists of propagating and evanescent waves. Information from details which are smaller than half a wavelength is contained in the evanescent fieldl. A classical optical microscope detects only propagating waves and this is why lateral resolution is limited to about half wavelength of the light source. Near field optical technique' is based on detection of propagating and evanescent waves in the near field of the sample i.e. at a distance far less than the wavelength À.

The probe detects the light by collecting a part of the propagative wave and by frustrating a part of the evanescent field laying near the surface of the sample. In order to reach a good lateral resolution the apex of the probe must have a radius of curvature smaller than 100 nm. The probe is here a multimode optical fiber chemically etched by HF (Hydrogen Fluoride). This kind of microscopy is known as Scanning Near Field Optical Microscopy (SNOM). Our experimental set-up that we had also described in ref. [START_REF] Berguiga | [END_REF].

With this microscope we are able to detect two different signals which are respectively optical signal and interaction force between probe and sample surface well known as shear force interaction. Since each signal could be chosen as the signal of feedback loop, the microscope has two working modes. But in the following section we will present only one mode whose shear force interaction is the signal of feedback 100p4.5.

The principle is quite similar to tapping mode in AFM, but here the probe vibrates parallel to the surface in one of its eigenmodes. When the tip approaches the surface, tip oscillation is damped due to lateral atomic force known as shear forcés

The effect of this force starts at ten nanometers above the sample surface (actually this starting value oscillates between 5 to 20nm), and increases continuously when distance between surface and sample diminishes.

The feedback loop locks the measured force in a fixed value, which corresponds to a fixed distance between sample surface and tip. And so when the tip scans an area of the surface, vertical motion of tip to maintain the distance constant gives the sample topography. The vertical resolution reached is below a nanometer and lateral resolution is in the dimension order of the tip end.

As the tip is very near from the surface, (we can set the distance from 1nm to IOnm), simultaneously the tip is able to detect the evanescent field emitted by the surface. The microscope turns out simultaneously two images: a topographic and an optical image.

In order to generate evanescent field, a 633 nm He-Ne laser illuminates a semi-cylindrical prism in total internal reflection configuration. The aged fiber was glued on a coverglass with transparent glue. A gel ensures optical contact between the coverglass and the prism. So evanescent field and propagating field are partly collected by the fiber probe and transmitted to a photomultiplier-tube. After that the optical signal is amplified and acquired by a Nanoscope E Digital Instrument system.

Low coherence reflectometry

The low coherence reflectometry technique is based on white light interferometry. Optical wave propagation defects are sligthly reflective and can be detected along single mode waveguides. We designed our system to derive both the amplitudes and the phases of the interferograms. We can analyze a fiber sample length of 2 meters, detect and localize two reflecting defects with a high accuracy and measure propagation and reflection parameters. We applied this technique and used our system on G652 fibers. The optical system is a Michelson single mode fiber-interferometer aligned around a 2-by -2 coupler. The broadband infrared light source ( 1.3 or 1.5 urn incoherent light emitting diode or amplified spontaneous emission fiber source) and the photodetector are respectively connected to the input and ouput arms of the coupler. The sample under test is connected to the test arm fiber. The wave propagating along the reference arm is collimated at the output of the reference arm fiber on a parabolic mirror, reflected-and-folded on a comer cube and reflected on a fixed plane mirror back to the coupler. The reference arm comer cube stands on an airborne translation stage back-to-back with two comer cubes. These two reflectors are the part of the variable length arm of a second Michelson interferometer associated to a frequency stabilized HeNe laser fringe counter. Under these conditions the HeNe fringe detector counts one fringe every 1/8 of a HeNe wavelength, and the optical path lengths of both red and infrared reference arms are known with a very high accuracy. The red fringe signals reshaped through the counter are used to sample the infrared interferograms which are processed afterwards.

Four to five fiber samples ( 50 ± 0.1 em) of each fiber under test have been connected to our measuring system. Samples of fibers which were not put under stress were used as references. The minimum detectable power of our system is approximately 5 pW when we use a stable and high power light source (20 mW). Its dynamic range is 75 dB, which means that a 6.10-4 difference between two refracted indices across a plane can be detected. We process the signals reflected by the cleaved output faces of the fiber samples. Each fiber sample is probed ten times: for each sample ten interferograms are detected and stored. The current injected in the source is kept constant and its light spectrum is assumed to remain the same. The interferograms Fourier transforms are averaged: both the amplitude and the phase are derived. The amplitude depends on many parameters-connectorsinsertion losses, fiber output angle and pollution which cannot be controlled. The phase then seems to be a far more reliable and appropriate parameter than the amplitude.

Mechanical measurements

Measurementshave been implemented according to the international standard lEe 793 12. Tested fiber length is one meter.

The static fatigue parameters were measured by a static bending test in which fiber is winded around a calibrated ceramic mandrel. Fibers are dried at room temperature and ambient atmosphere before testing. Fiber failure is measured optically. Mandrel diameters were 2.7, 2.8, 3.0 and 3.2 mm, which corresponds to stress values of 3.34, 3.22, 3.0 and 2.8GPa respectively.

PREPARATION OF FIBER FOR OBSERVATION

Aging conditions

We have chosen commercial optical fiber G652. They were immersed in deionized water at 65 and 85°C. Time aging are three, six, nine months and one year. Fiber was furnished and aged by C.N.E.T Lannion (France). The process aging at 85°C should be different than at 65°C because at 85°C the integrity of the coating is altered. We have noticed the color change of polymer at 85°C for three, six and nine months aged fiber.

Stripping technique

If we want to study silica surface of aged optical fiber, the polymer coating should be removed without degrading the mechanical properties of fiber. There are three ways to strip coating from the fiber. The most used technique is mechanical stripping with strippers, but it damages the fiber surface. The second way, which is the easiest, is to swell the coating using a solvent such as methylene chloride. Since the coating loses its adherence on the cladding, it becomes very easy to remove it. In the last solution, the polymer coating is stripped by an acid attack", The acid is usually sulfuric acid at high temperature.

We have chosen acid stripping for several reasons. As we discussed previously, Rondinella and Matthewson7 studied the surface roughness of aged fiber and showed it did not reveal the aging when the fiber is stripped by methylene chloride. In contrast with acid stripping, the surface is correlated with the residual strength. They interpret this fact by the existence of a hydrated surface layer of silica gel that would be remove only by acid. While the fiber is aged, the penetrating water through coating to the silica surface may produce this layer. Moreover it was demonstratedin reference 6) coating was stripped in hot concentrated sulfuric acid, the strength ofthe fiber is not affected.

In a former3 work, the stripping conditions (see table I) were not as much efficient as we expected to eliminate totally the polymer cladding. Despite the increase of the immersion time in acid, elements of polymer still remains with an average lateral size of 3 to 5 hundreds nanometersand one hundred nanometer high.

This observation led us to prepare harsher stripping process. Two different processes (Cf table I) were set-up. These two processes consist to decompose stripping in two step. The first one suppresses the majority of polymer cladding by acid attack for process A or by swelling with a solvent for process B. In the second step, hot acid (l80-200°C) should attack polymer which has resisted to the first step.

A study of the stripping of each process has been done versus the time immersion in acid (lmn to 20mn). From images obtained with SEM (Cf Figure I) and with the Shear Force Microscope (SFM) in figure 2. We can remark that:

-Even if process A and B are harsher, polymer residues still remains.

-The increase in time of immersion in acid does not improve at all the efficiency of stripping. In other words size and density of residues on the surface is quite similar.

If we compare these processes with the process of former work', process A and B are better. Per unit of surface, residues are less numerous than with former process, and the average lateral size and height of residues decrease respectively to 100 nm and 20 nm. This difficulty to eliminate totally the polymer coating could be explained by a strong adherence of the polymer coating on the silica surface of the present optical fiber.

Former Process Immersionin H2S04 at 160°C+ rinsing in deionized water 30s

Process A Immersion in H2S04 at 180-20QoC+Swelling with methylene chloride + rinsing in deionized water 30s +Immersion inH2S04 at 180-20QoC+rinsing in deionized water 30s

Process B

Swelling with methylene chloride+ rinsing in deionized water 30s+

Immersion in H2S04 at 180-200°C+rinsing in deionized water 30s For relatively important domain of (lO!-lm)2, we observe still polymer residues, but domains of the order of (1 ~/ are often exempt of residues. In the next section we will present images and statistical results of small clean area of optical fibers with aging conditions described in paragraph3.1 and stripped with process B.

RESULTS AND DISCUSSION

Study of topography of aged fiber

The surface topography of aged fiber at 65°C and 85°C during 3, 6, 9, 12 months has been observed with our near field microscope. Some typical images of these surfaces are presented in figures 4 to 9. No apparent modification of the surface with aging in harsh environment was observed. Moreover it is rather hard to interpret these images only with qualitative observations. Table II shows the roughness RMS versus aging conditions and it confirmsthat none evolution is noticeable.

Since the number of images of fiber surface for each condition of aging is important, statistical study is rather tedious with only the roughness parameter. We use the Power Spectral Density (PSD), which is the frequency spectrum of the surface roughness measured in inverse length units". The PSD provides an information about roughness at different scales of a topographic image. For example, ifthe surface presents a periodic structure or a change in the structure of the surface versus or spatial frequency, this behavior should be detected with the PSD graph.

For each condition of aging, PSD of many scanned surfaces of I to (2 !-lm2) are compared. We note that for fiber aged 3 months at 85°C, PSD curves are very narrow to each other (see figure 10).This point indicates that structure of fiber surface is homogenous. Same feature for fiber aged 3, 6, 9 months at 65°C and 85°C was observed. About fiber aged 12 months at 85°C (figure 11), a dispersion in the PSD graphs occurs at high spatial frequency. Still in figure 14, this dispersion raises when we go to higher spatial frequency.Fibers aged 12 months at 65°C seem to be an intermediate case where the apparition of dispersion (Cf. fig. 12) is less important than in figure 11. In the graph of figure 11 (so for fiber aged 12 months at 85°C) the starting point of dispersion is for lower spatial frequency than in figure 12 (i.e. fiber aged 12 months at 65°C). We can interpret these results by two phenomena:The surface begins to change not at every scale but for the smallest details. At this scale the roughness increases but heterogeneously, in other words some areas have modified structure for fine details whereas other areas are less affected. When this effect becomes important for the fine details, the change in the structure of the surface for greater details starts and involve heterogeneously at this scale too. And so on for greater scale... These results show despite harsh aging conditions, the evolution of fiber surface is very weak and surface is very similar to non aged optical fiber and thin differences are observable for one year aging and at nanometric scale. Those observations are not in contradiction with mechanicaltests and low coherence reflectometry.

Optical information

Near field optical images which are obtained do not reveal any particular information about aging of optical fiber. Optical image could give a complementaryinformation about surface and subsurface of the fiber that topography is not able to detect. To illustrate this affirmation, on a surface of a fiber aged at 9 months., a localized structure along a line of few micrometers length and 200 nm wide went across the optical image and did not appear on the corresponding topographical image. It turns out that this detected structure can be interpreted as a local refractive index variation or a failure in the silica subtract which is not resolved by the SFM. When we compare this optical image and numerical simulationll, the nearest result is for a crack of 20 nm deep and 10 nm wide.

Low coherence reflectometry

Measurements have been implemented according to the conditions described above. All samples have been probed along their total lengths and no characteristic echo has been detected. In other words, the aging conditions considered did not induce defects with such abrupt variations of the refractive index.

The phase of the signal is a function of the wave number and we calculate its derivative. Two sets of derivatives are plotted against the wave number on figure 13. They correspond to two different types of fibers and two aging test durations at 85 "C. These curves have a parabolic shape and present an extremum at the fiber zero dispersion wave number (ZDWN). We will focus on this last parameter which seems to depend on the fiber intrinsic characteristics only. We experimentally observe a variation of 10') mu" this zero dispersion wave number for different samples of the same fiber which we consider as the lowest detectable limit with our measurement technique. This sensitivity is high enough to detect a noticeable difference on the evolutions of the two types of fibers. More precisely the ZDWN of the fiber with the «standard»protection coating is approximately 0.758 mu' and the ZDWN of the fluoride protection coating is close to 0.752 mu'. These two values correspond to two different fiber structures as we could observe in SEM pictures. This demonstrates that we can quantitatively characterize fibers and identify them without ambiguity Aging tests at 85 "C have been applied to the fibers with the «standard» coating. The phase derivatives vs wave number are plotted on figure 13b for several test durations. We see no difference between the curves and this feeling is confirmed by the measurement of the ZDWN : they are in the order of 10-4 mu' which is not significant. We can draw the same conclusions from the results obtained on the fibers with the fluoride coating.

Static fatigue tests

The static fatigue measurements show that the failure time increases when the fiber is aged in hot water. This is exemplified by figures 14 to 17. This effect depends on the nature of the coating and also on temperature. The evolution of the failure time is not monotonous. It is maximum after 3 to 6 months, then it decreases beyond the initial value. Note that the behavior at 65 "C is different with the fluorinated coating:. In a first step, failure time decreases, then it is enlarged. Coating seems damaged after 15 months at 85 "C.

These results may appear paradoxical as water is known to decrease fiber strength and reliability. The overall process is complex as it involves interactions between water molecules and silica, but also polymeric coating.. Among other possibilities, water could reinforce the bonding adhesion between coating and glass surface.

CONCLUSION

Silica fibers have been aged in water at 65 and 85 oc. No significant change was observed at fiber surface after special stripping. No aging induced stresses was detected by low coherence reflectometry. Static fatigue tests shows that failure time is larger for aged fibers after 3 to 6 months. For larger aging times, fiber strength becomes much smaller, which may be correlated to the degradation of the coating. ...

2.0

X 0.500 PM/div 2 10.000 """diu

Figurel:Figure 3 :Figure 4 :Figure 6 :

 346 Figurel: SEM image of an optical liber stripped withProcess A,

Figure 8 :Figure 5 -Figure 7 :Figure 9 :figure 10 :figure 11 :figure 12 :

 8579101112 Figure 8: Image ofthe surface of 12 months aged at 65°C

Figure 17 :

 17 Figure 17 : Fiber 2 aged at 85°C. Right scale correspond to the 3.2 mandrel.

Table I :

 I Process of stripping

Table II :

 II RMS roughness for IJ.1m2area versus aging conditions of optical fiber.

		Non aaed fiber	3 months	6 months	9months	12 months
	16Soc 18Soc	0.5nm	0.5nm OA-0.5nm	OAnm O.3-0.4nm	0.5nm No measurement	OAnm 0.3-0.8nm


Nathalie Huet1

Eric Tanguy2
email: eric.tanguy@physique.univ-nantes.fr

Armand Vinçotte1

Vassilis Zafiropulos3

Using lasers for cleaning ceramic and plaster patrimonial objects

Traditional cleaning techniques (sandblasting, chemical, ... ) turn out to be ineffective or damaging on some ceramics and plaster objects. It's especially true on outdoor ceraIJÙCScovered with black crust or indoor biscuits wares and plasters exposed to dust and repeated handling. The use of lasers to remove dirt, such as for stone, has then been considered. This paper study the behaviour of these materials under the impact of different type of lasers irradiation (Nd:YAG first and third harmonic, excimer, ... ). According to the type of dirt and underlying material, the effectiveness of cleaning depends on the laser wavelength and its energy density. The results are satisfactory with no damage and no coloration change on plaster cleaning by Nd: yAG third harmonic.

INTRODUCTION

The use of Nd:YAG lasers to remove dirt from stone monuments is gradually replacing [il the more abrasive technique of sandblasting. This trend has caused Arc' Antique to consider lasers as a solution for the cleaning of outdoor glazed ceramics (stoneware bowls) covered with a blackish layer of pollution (Figure 1). In fact, the standard techniques using chemical products have proved totally ineffective on these objects and sandblasting is relatively damaging to the glaze. The promising results obtained when cleaning glazed ceraIJÙCS [ii] suggested to us that the technique could be also used for some indoor ceramics, like biscuit ware, and plasters. These materials are indeed often covered in museum dirt (dust, grease, etc.) that is difficult to remove without damaging the object.

THE CHOICE OF THE WAVELENGTH

Following initial tests with an infrared Nd:YAG laser, it was decided that research have to be undertaken into the behaviour of ceramic and plaster materials under the impact of lasers, into the types of laser most suited to the removal of the most commonly found dirt, and into optimisation of the different operational parameters. The first study started in FORTH (Crete), where we were able to take advantage of the expertise of the ARTT team in the cleaning field and to use their high performance equipment (different lasers, on-site monitoring techniques). It consisted of identifying the materials cleaned and layers of dirt to be removed. We had the opportunity to tryout the available cleaning options for different categories of dirt (pollution, paint, dust, oxide staining, grease) on several types of material (pottery, porcelain, biscuit ware, plaster) under the impact of different types of laser (Nd:YAG first and third harmonic, excimer). It was demonstrated that the effectiveness of this cleaning is complex, as, for each type of dirt and underlying material, the behaviour differs depending on the type of laser used (and therefore the wavelength of the laser) due to the physical mechanism of ablation involved (photothermal/photochemical).

Some results appeared really encouraging, such as the removal, with an excimer laser (KrF, À = 248 nm), of layers of dust from "rough" porous and fragile objects, such as biscuit or plaster, with no damage. Similarly, spectacular results can be obtained with a third harmonic Nd:Y AG (À = 355 nm) in removing accumulated dirt caused by external pollution from certain glazed ceramics (Figure 2) or dust on rough clay mass and plaster. Highly enlarged images demonstrate that the surface is not damaged, unlike with standard techniques such as micro sandblasting.

PROBLEMS AND CAUSES

Despite these good results, certain tests were catastrophic. A number of observations are leading us to make a more detailed study of phenomena using SEM, EDS, X-Ray diffraction, infrared spectrometry. The first problem is due to the thermal effect during the use of Nd:YAG lasers (À = 1.06 11m)provoking cracking and fusion of glazes that can be avoided with a precise control of the fluence. But the major effect is a change of coloration which affect glazes, painting, clay mass and plaster. Looking closely at these disappointing tests, we tried to analyse the causes or at less the material on wich these phenomena appeared. Two problems were point out on certain materials depending of the wavelenghts used: a yellowing [iii] and colour alterations. With the infrared Nd:YAG laser, usually used for stone, the surface becomes yellowish (Figure 2) or greyish. If on stone it could be identified as patina, the same explanation is not available on ceraIJÙcs where ageing cause no yellowing. It seems that this colour is not caused by an alteration of the material but rather a layer of material (dirt") redeposited or still deposited on the surface. Although, as many others researchers, it has not yet been possible to identify this layer even using sophisticated techniques (ESCA). Nevertheless, a second cleaning with another type of laser, such as the third harmonic Nd:YAG, can totally remove it and restore the


Ke Wu

Huazhong Shu

Jean-Louis Dillenseger
email: jean-louis.dillenseger@univ-rennes1.fr

Region and boundary feature estimation on ultrasound images using moment invariants

Keywords: Texture analysis, ultrasound images, speckle, moment invariants, min-cut/maw-flow graph cut

Region and boundary feature estimation on ultrasound images using moment invariants.

Introduction

Segmentation of ultrasound images (US) is a difficult task because their characteristic feature is speckle. Speckle arises from signal interference caused by tissue micro-inhomogeneities (tissue cells, capillaries, blood cells, etc). This coherent summation of backscattered signals forms a spatial distribution of speckle that is specific to the density and distribution of the scatterers and thus to the nature of the tissue [START_REF] Wagner | Statistics of speckle in ultrasound B-scans[END_REF]. Several speckle spatial distribution models have been proposed in the literature according to the level of organization and density of scatterers within the tissues (Rayleigh-, Rician-, K-, Nakagami-distributions, etc.). Hence, speckle has been used in several ways to characterize tissues [START_REF] Noble | Ultrasound image segmentation: a survey[END_REF]. Some authors have tried to estimate the speckle distribution model using statistics on gray level intensities (e.g. [START_REF] Seabra | Rayleigh mixture model for plaque characterization in intravascular ultrasound[END_REF]), but this class of methods overlooks the spatial nature of the distribution. Texture analysis has been traditionally used to characterize the spatial distribution of patterns. Besides studies that aimed to directly estimate the speckle distribution (see [START_REF] Shankar | The use of the compound probability density function in ultrasonic tissue characterization[END_REF] for example), the use of general texture analysis methods had some success in ultrasound image segmentation (see [START_REF] Noble | Ultrasound image segmentation: a survey[END_REF] for a survey). However, texture analysis needs relatively large windows to perform feature estimation; this leads to a lack of precision, especially at tissue boundaries. Other US image characteristics can also make feature extraction more problematic. Because of the ultrasound beam form, the size (or scale) of the speckle increases according to the distance from the ultrasound probe. Moreover, for circular probes, the ultrasound beam directions radiate from the probe center because of its geometry. These radial directions have a direct impact on the speckle orientations, which are different depending on the speckle position across the image (see the circular speckle organization in Fig. 1). In these cases, a texture analysis tool should be able to extract features that are invariant with regard to both scale and orientation. Some attempts have been proposed to extract invariant texture features using Gabor filters [START_REF] Zhang | Invariant texture segmentation via circular Gabor filters[END_REF][START_REF] Han | Rotation-invariant and scale-invariant Gabor features for texture image retrieval[END_REF]. In ultrasound image analysis, multiscale and multiorientation twodimensional (2-D) Gabor filter banks (composed of six different orientations and two different scales) have been used to characterize the prostate boundaries and incorporate into a Kernel Support Vector machine for texture differentiation [START_REF] Shen | Segmentation of prostate boundaries from ultrasound images using statistical shape model[END_REF][START_REF] Zhan | Deformable segmentation of 3-D ultrasound prostate images using statistical texture matching method[END_REF]. The likelihood obtained by this means is then used to drive a deformable surface model; however, this scheme is relatively complex.

Moment-based texture analysis has been introduced in the literature. The main idea of this method is the extraction of some features by calculating a set of moments based on regions of interests (ROIs) on texture images. For example, low order geometric moments correspond to geometric meaning of texture in the ROI, such as mass, centroid and principal axes. Other classes of moments are not directly related to any specific properties of an image, but their suitability for texture characterization has been confirmed [START_REF] Bigun | N-folded symmetries by complex moments in Gabor space and their application to unsupervised texture segmentation[END_REF][START_REF] Tuceryan | Moment-based texture segmentation[END_REF][START_REF] Mandal | Image indexing using moments and wavelets[END_REF]. In his pioneer work, Hu proposed to derive a set of seven 2-D moment invariants from the classical geometrical moments according to translation, rotation and scale [START_REF] Hu | Visual pattern recognition by moment invariants[END_REF]. Hu's moments have been widely used for pattern recognition, but could also be used for texture analysis.

Moments with orthogonal basis functions (e.g. Zernike, pseudo-Zernike polynomials) can represent an image by a set of mutually independent descriptors, and thus have a minimal amount of information redundancy. They were shown to have better image description capability and to be more robust to noise than geometric moments. In the past few years, the construction of orthogonal moment invariants with respect to geometric transformations such as translation, scale and rotation has been investigated [START_REF] Chen | Combined invariants to similarity transformation and to blur using orthogonal Zernike moments[END_REF][START_REF] Zhang | Object recognition by a complete set of pseudo-Zernike moment invariants[END_REF][START_REF] Zhang | Construction of a complete set of orthogonal Fourier-Mellin moment invariants for pattern recognition applications[END_REF]. These moment invariants have been applied to pattern recognition. In this paper we propose to use them for texture analysis and aim to identify some of their properties that can be used in a simple segmentation scheme. The organization of this paper is as follows. Section 2 presents the three different sets of moment invariants. The performance of the features extracted by the moment invariants on simulated and real data is then evaluated in section 3. This same section presents also some gional or boundary information. Several segmentation methods applied on gray scale images use together regional and boundary information (active contours, graph-cut, etc.). With the moment invariants, these techniques could be applied on US images. At the end of section 3 we will try to prove this assessment on a min-cut/max-flow graph cut scheme.

The originality of this paper is to highlight the usability of moment invariants to extract texture features when the texture is varying in scale and orientation as in US images. Some authors used Gabor filters bank [START_REF] Zhang | Invariant texture segmentation via circular Gabor filters[END_REF][START_REF] Han | Rotation-invariant and scale-invariant Gabor features for texture image retrieval[END_REF] After some test we were surprised to see the property that the moments of repetition q=0 provide regional information and these of repetition q=1 provide boundary information. This fact seems for us sufficiently interesting for a publication because several methods applied on gray scale images use regional and boundary information (active contours, etc.). With the orthogonal moment invariants, these techniques can now be applied on US images. We tried to prove this assessment on the graph cut scheme, but the graph cut itself is not here the main purpose of the paper.

Moment invariants

As specified in the introduction, the texture analysis on an ultrasound image should be performed using some descriptors that are invariant to rotation because of the radial organization of the speckle and invariant to image-scaling because the speckle size is directly related to the ultrasound beam diameter which is not constant over the image. In this study, we used complete sets of moment invariants reported in [START_REF] Chen | Combined invariants to similarity transformation and to blur using orthogonal Zernike moments[END_REF][START_REF] Zhang | Object recognition by a complete set of pseudo-Zernike moment invariants[END_REF][START_REF] Zhang | Construction of a complete set of orthogonal Fourier-Mellin moment invariants for pattern recognition applications[END_REF]. The following subsections give brief descriptions of these invariants.

Zernike moment invariants

The Zernike moment formulation is one of the most popular orthogonal moments. In the polar coordinate system, the 2-D Zernike moment (ZM ) of order p and repetition q (where p ≥ 0, |q| ≤ p and p -|q| = even) of an image intensity function f (r, θ) is defined as in [START_REF] Mukundan | Moment functions in image analysis: theory and applications[END_REF]:

ZM f p,q = p + 1 π ∫ 2π 0 ∫ 1 0 R p,q (r) exp(-jqθ)f (r, θ)rdrdθ (1)
where R p,q (r) is the orthogonal radial polynomial:

R p,q (r) = p-|q| 2 ∑ k=0 (-1) k (p -k)! k! ( p+|q| 2 -k ) ! ( p-|q| 2 -k ) ! r p-2k (2)
When an image undergoes a rotation of angle α, the Zernike moments of the transformed image are given by ZM f p,q exp(-jqα). Thus, the magnitude of the Zernike moment, i.e. |ZM f p,q |, is invariant to rotation. However, as indicated by Flusser [START_REF] Flusser | On the independence of rotation moment invariants[END_REF], the magnitudes do not generate a complete set of invariants, and completeness is an important property for assessment of image descriptors. Additionally, the original Zernike moments are not invariant to image scaling.

To solve this problem, we proposed an approach to derive a complete set of Zernike moment invariants (ZM I) with respect to rotation and scaling of order p and repetition q (where p = q + 2m with m ≥ 0) [START_REF] Chen | Combined invariants to similarity transformation and to blur using orthogonal Zernike moments[END_REF]:

ZM I f p=q+2m,q = exp -jqθ f m ∑ k=0 m ∑ l=k Γ -(q+2l+2) f c q m,l d q l,k ZM f q+2k,q (3)
where ZM f q+2k,q is the original Zernike moment defined in (1), θ f = arg(ZM f 1,1 ), Γ f = √ ZM f 0,0 and c q m,l and d q l,k are given by: c q m,l = (-1) m-l q + 2m + 1 π

(q + m + l)! k!(m -l)!(q + l)! d q l,k = π l!(q + l)! (l -k)!(q + l + k + 1)!
Equation [START_REF] Seabra | Rayleigh mixture model for plaque characterization in intravascular ultrasound[END_REF] shows that the moment invariants are a linear combination of their original moments, thus they maintain the capacity for texture characterization.

Pseudo-Zernike moment invariants

The 2-D Pseudo-Zernike moment (P ZM ) of order p and repetition q (where |q| ≤ p) was then defined as in [START_REF] Teh | On image analysis by the methods of moments[END_REF]:

P ZM f p,q = p + 1 π ∫ 2π 0 ∫ 1 0 P p,q (r) exp(-jqθ)f (r, θ)rdrdθ (4)
where P p,q (r) is the orthogonal radial polynomial given by:

P p,q (r) = p-|q| ∑ k=0 (-1) k (2p + 1 -k)! k! (p + |q| -k + 1)! (p -|q| -k)! r p-k (5)
Pseudo-Zernike moments have been proven to be superior to Zernike moments in terms of their feature representation capabilities and robustness in the presence of image quantization error and noise [START_REF] Teh | On image analysis by the methods of moments[END_REF].

Following the same methodology as for the Zernike moments, a complete set of pseudo-Zernike moment invariants (P ZM I) with respect to rotation and scaling of order p and repetition q (where p = q + m with m ≥ 0) was constructed as in [START_REF] Zhang | Object recognition by a complete set of pseudo-Zernike moment invariants[END_REF]:

P ZM I f p=q+m,q = exp -jqθ f m ∑ k=0 m ∑ l=k Γ -(q+l+2) f c q m,l d q l,k P ZM f q+k,q (6)
where P ZM f q+k,q is the original pseudo-Zernike moment,

θ f = arg(P ZM f 1,1 ), Γ f = √ P ZM f 0,0
and c q m,l and d q l,k are given by:

c q m,l = (-1) m-l q + m + 1 π (2q + m + 1 + l)! l!(m -l)!(2q + 1 + l)! (7)
d q l,k = 2π l!(2q + l + 1)! (l -k)!(2q + l + k + 2)! (8)

Orthogonal Fourier-Mellin moment invariants

The 2-D orthogonal Fourier-Mellin moment (OF M M ) of order p with repetition q of an image intensity function f (r, θ) was defined as in [START_REF] Sheng | Orthogonal Fourier-Mellin moments for invariant pattern recognition[END_REF]:

OF M M f p,q = p + 1 π ∫ 2π 0 ∫ 1 0 Q p (r) exp(-jqθ)f (r, θ)rdr (9)
where Q p (r) is a radial polynomial given by:

Q p (r) = p ∑ k=0 (-1) p+k c p,k r k (10) c p,k = (-1) p+k (p + k + 1)! (p -k)!k!(k + 1)! (11)
For region of small radial distance the description by orthogonal Fourier-Mellin moments has been found to be better than that by the Zernike moments in terms of image-reconstruction errors and signal-to-noise ratio [START_REF] Sheng | Orthogonal Fourier-Mellin moments for invariant pattern recognition[END_REF]. It can be proven that when q = 0, OF M M f p,0 = P ZM f p,0 .. A complete set of orthogonal Fourier-Mellin moment invariants (OF M M I) with respect to rotation and scaling of order p and repetition q of an image intensity function f (r, θ) was constructed as in [START_REF] Zhang | Construction of a complete set of orthogonal Fourier-Mellin moment invariants for pattern recognition applications[END_REF]:

OF M M I f p,q = exp -jqθ f p ∑ k=0 p + 1 k + 1 ( p ∑ l=k Γ -(l+2) f c p,l d l,k ) OF M M f k,q (12)
where OF M M f k,q is the original orthogonal Fourier-Mellin moment,

θ f = arg(OF M M f 1,1 ), Γ f = √ OF M M f 0,0 , c p,l
is defined by ( 11) and d l,k is given by:

d l,k = (2k + 2) l!(l + 1)! (l -k)!(l + k + 2)! (13)
As for the original moments, it can be proven that when q = 0, OF M M I f p,0 = P ZM I f p,0 .

Results and discussion

In this section we evaluate the behavior of the three moment-invariant sets relative to their original sets, as well as to the Hu's moment invariants [START_REF] Hu | Visual pattern recognition by moment invariants[END_REF]. Because the magnitude of the orthogonal complex moments is supposedly invariant to rotation, we will only use this modulus in the rest of the paper. Hu's moments will be denoted as Hu n with n indicating the same invariant moment number as given in Hu's paper.

Material

The evaluations were performed on a synthetic image and on in vivo ultrasound data.

The synthetic image was composed of an elliptically-shaped region set within a background region. The ultrasound simulation was performed by the method described in [START_REF] Dillenseger | Fast simulation of ultrasound images from a CT volume[END_REF]. This method predicts the appearance and properties of a B-Scan ultrasound image from the spatial distribution of point scatterers. In our case, each region was characterized by a specific acoustic impedance and spatial distribution of scatterers. The parameters defined in [START_REF] Dillenseger | Fast simulation of ultrasound images from a CT volume[END_REF] for fat and liver were used to differentiate the two regions. This scatterer model served as input for the simulation of an ultrasound image acquired with a C2 5 circular ultrasound scanning 3.5 MHz probe (Fig. 1-a). In this image the gray scale depended only on the spatial distribution of the scatterers. The image size was 512 × 512.

The in-vivo data was a transrectal prostate image acquired intraoperatively on an Ablatherm device (Fig. 1-b). The transrectal ultrasound imaging probe operated at 7.5 MHz. A slice had 500 × 490 pixels with a transverse pixel size of 0.154 mm/pixel and a thickness of 2 mm.

Invariance to rotation and scale

First, we evaluated the rotation and scale invariance of our moment-invariant sets. We extracted a 51 × 51 ROI window from the bottom of the simulated ultrasound image (see Fig. 1-a) and applied rotation and scaling transformations to this sub-image. The moment sets were then computed in the center of this transformed sub-image.

Rotation invariance: the extracted sub-image was first rotated by an angle from 0 to 90.

We evaluated the variation of the modulus of the moment sets and the moment invariants with regard to the rotation angle. In this paper, we present this variation and compare the sets of classical moments with their invariant forms for P ZM and P ZM I only. However, these results can be generalized to all the other moment sets. So, as an example, Fig. 2-a shows the variation of the moduli of P ZM 3,0 and P ZM I 3,0 in relation to rotation angle. Surprisingly, even if the moduli of the classical moments (ZM , P ZM and OF M M ) should be theoretically invariant to rotation they show some relatively high variations around a mean value when the angle changes.

Our proposed invariant sets (ZM I, P ZM I and OF M M I) remains almost constant when the angle changes.

Scale invariance: we tested this invariance by increasing the moment window size from 7 × 7 to 31 × 31 pixels. The same speckle pattern was seen at several scales when the window size was

varied. An example of moment behavior with scale can be seen in Fig. 2-b. In this figure, as was also the case for all the moment sets, the value of the classical moments is window-dependent but remains stable for the moment invariants.

Moment widow size

Although the moment invariants are scale independent, window size is an important parameter for any texture feature extraction method. A compromise must be found between a large window size to capture the feature and a small window size to accurately locate the information. Window size is usually related to the spatial scale of the texture. We choose to illustrate the effect of the window size on the estimation of the P ZM I 1,1 feature of the synthetic image. This feature presents some fine structures which will directly impacted by the window size. In this image, the speckle longitudinal size was about 20 pixels. Fig. 3 shows the P ZM I 1,1 feature computed with a window sizes of 11 × 11, 21 × 21 and 31 × 31. It can be clearly notified that the feature homogeneity increases with the window size but at the expense of the thinness of the features.

In this case, the 21 × 21 window seemed to be a good compromise between the homogeneity of the features and the accuracy to extract some image components. This window size will be used in the rest of this paper. This behavior was confirmed for the other moment feature measurements.

Order and repetition characteristics

In this section we will try to evaluate the behavior of the moment invariants according to the order p and repetition q parameters.

Order itself has a relatively low influence on the final result. As an example Fig. 4 shows the resulting images for P ZM I p,0 with p increasing from 1 to 3. We can observe that the images are relatively similar. Maybe the images seems to be a bit more irregular for higher orders but the details like then contours are more well defined.

The effect of the repetition is more interesting in the context of further segmentation. The behavior of the repetition can be clearly seen on Fig. 5. We can note that Hu 1 and the moments of repetition q = 0 (see first column of Fig. 5) are able to extract some features that are sensitive to regional texture information. On all these images, it can be seen that the moments and the moment invariants are able to characterize differently the two textured regions. However, when q = 1, and also to a certain extent for Hu 2 (see second column of Fig. 5), the features are sensitive to rupture between two textures. These later features seem to be very informative about the boundaries of the objects and could be similar to a "gradient of texture".

This regional and boundaries extraction behavior can also be seen on the features estimated from the real transrectal prostate ultrasound data (see Fig. 6 for the pseudo-Zernike case). On this figure, the features extracted using the moment invariants are also more homogeneous than those obtained using the classical moments.

Invariance properties for texture estimation

We tried to prove the invariance properties of moment invariants for texture feature estimation both qualitatively and quantitatively. We made these tests on the simulated ultrasound image in Fig. 1-a. Qualitative results can be seen on Fig. 5 which shows a selection of extracted features from several kinds of moments and moment invariants with different values of order p and repetition q. For the quantitative evaluation, we choose forty points (the crosses on Fig. 1a) in the same region but at different locations, thus having different speckle orientation and scale. The measurement of the Relative Standard Deviation (RSD -the absolute value of the ratio of the standard deviation to the mean) of the features computed on those points was used to estimate the invariance ability (Table 1).

It can be seen on Table 1 that for the classical moments the value of the texture feature is very dependent on the location of the ROI because of differences in speckle orientation and scale. In contrast, moment invariants showed a relatively (or more) stable measure. The proposed orthogonal moment invariants also showed more stable features than Hu's moment invariants. This behavior is also noticeable on the feature images of Fig. 5 when visually comparing the images of the classical moments and its equivalent (same order and repetition) moment invariants. In all the cases the features in a same region are more homogeneous when estimated by moment invariants. This is specially true for P ZM and OF F M which seems to have an incoherent behavior in a region sharing the same speckle distribution. This could be generalized for moments of any order or repetition.

However, it could be noted that the RDS for the moments (classical and invariants) with repetition q = 1 and for Hu 2 were much higher than for the moments with q = 0 and for Hu 1 . This increase is explained by the fact that in these cases where the moments features are sensitive to the rupture between two textures, the means of the values were close to 0.

In conclusion, the features extracted by moments invariants, are relatively homogeneous (small RDS) within a region sharing the same speckle distribution. They can so be used as texture descriptors. In all the cases, moments invariants showed a more homogeneous behavior than the classical moments.

Proposal for using moment texture features

In texture analysis, a set of texture features (moments and others) are usually measured on the image, and a classification process is then performed in order to characterize the objects of interest. However, as the moment-invariant features with repetition q = 0 provide regional information and those with q = 1 provide information on the boundaries of a textured region, this allows us to use the texture features differently. In many segmentation techniques, boundary information alone or in combination with regional information is needed, i.e. active contours, deformable models, level sets, etc. Some of these techniques have been adapted to deal with ultrasound images and require boundary information (e.g. [START_REF] Chen | An early vision-based snake model for ultrasound image segmentation[END_REF] where some Gabor based invariants where used to constraint snake models). Our two kinds of moment-invariant features (q = 0 regional information and q = 1 boundary information) could be directly used in such a segmentation scheme dealing with ultrasound images.

Graph-cut, a way to integrate region and boundary information

The aim of this section is to try to prove the previous assessment, without loss of generality, by using the two kinds of moment-invariant features in a min-cut/max-flow graph cut algorithm.

So we will show how the features extracted by orthogonal moments can be directly used in a combined region and boundary based segmentation algorithm and how the boundary information can improve the segmentation. However, the goal of this section is neither to evaluate finely the impact of the several parameters used in the graph-cut segmentation nor to compare its performance to other techniques.

The min-cut/max-flow graph cut algorithm popularized by Boykov is able to handle regional and boundary information simultaneously [START_REF] Boykov | Graph cuts and efficient N-D image segmentation[END_REF]. To summarize, the segmentation problem can then formulated by the energy function E T = E classif + λ • E continuity . E classif is an energy term coding the probability that a pixel belongs to the class "object" or "background" (a regional energy). E continuity is an energy term coding the degree of similarity or discontinuity between two neighboring pixels (a boundary energy). The coefficient λ controls the balance between these two energy terms. The energy terms are encoded within the graph as link weights estimated from the pixel values based on the region and continuity (or boundary) properties of the image [START_REF] Boykov | Graph cuts and efficient N-D image segmentation[END_REF][START_REF] Esneault | Graph cut liver segmentation for interstitial ultrasound therapy[END_REF].

To validate the usability of the moment invariants, we adapted a variant of the min-cut/maxflow graph cut algorithm presented in [START_REF] Esneault | Graph cut liver segmentation for interstitial ultrasound therapy[END_REF]. We used the moment-invariant features with repetition q = 0 in order to encode the regional information and the moment-invariant features with repetition q = 1 to encode the similarity between adjacent pixels. As an example, we will take the case with Zernike moment invariant features. In this method, the user labels interactively some pixels as "object" and others as "background" (e.g. Fig. 7-a or 8-a). The normalized histogram of these pixels in ZM I p,0 is used to define the probability density function P object (resp. P backg ). For a pixel u, the weights {u,S} (resp. {u,T})1 associated with the regional energy can be seen in Table 2. E continuity is encoded on the links between two adjacent pixels {u, v}. This weight should be high between two neighboring pixels within the same region and low on boundaries. We first normalized and inverted ZM I p,1 to obtain image ZM I ′ p,1 which was representative of the continuity between neighboring pixels. The weight between {u, v} was then given by:

ZM I ′ p,1 (u)+ZM I ′ p,1 (v) 2
(see Table 2).

We performed the evaluation for the Hu's moment invariants and the three orthogonal moment invariants on the simulated US image (Fig. 1-a). Some small areas of pixels were interactively selected in both the object and background regions (Fig. 7-a). These were used as seed points for the segmentation and also serve as training basis for the E classif -related initial weight assignment. The extracted region was compared with the real elliptically-shaped region.

Each pixel could be classified as true positive (with T P the number of true positives), false positive (resp. F P ), true negative (resp. T N ) or false negative (resp. F N ). The performance of our segmentation scheme is presented as fractional area difference

( F P -F N T P +F N )
, sensitivity

( T P T P +F N )
and accuracy

( T P -F P T P +F N
) .

The crucial points of our method are: the choice of the features used to encode the region and boundary information, and the balance coefficient λ.

After some exploratory work we noticed that the more accurate results were obtained by using relatively high order moments for the regional information Hu 1 , ZM I 2,0 , P ZM I 3,0 and OF M M I 3,0 and for the regional information Hu 2 , ZM I 1,1 , P ZM I 2,1 and OF M M I 2,1 . These features can be seen in Fig. 5 and are consistent with the remarks made in the previous section about more well defined details for higher moment orders.

The λ parameter setting was trickier because its influence is not linear and is data dependent.

First, in an exploratory study, we interactively tuned λ to analyze its influence. We noticed that, even with a relatively small λ, the continuity term was taken into account within the segmentation. Then, for a relatively large range of variation in λ, the segmentation results remained almost the same, with only some small changes on the boundaries. Finally, for very high values of λ (over 100) the segmentation accuracy decreased. Fig. 7 shows the impact of the variation of λ. The segmentation results presented were obtained using the ZM I set. For λ = 0, we only had the influence of the regional information of ZM I 2,0 (Fig. 7-b). With λ = 1, the boundary information of ZM I 1,1 was mixed with the regional information in order to delineate the region with a higher accuracy (Fig. 7-c).

Table 3 compares the performance between the segmentation results of the min-cut/maxflow graph cut algorithm with λ = 1, using methods based on Hu, Zernike, pseudo-Zernike and orthogonal Fourier-Mellin moment invariants and the equivalent classical moments. In all cases (except one), moment invariants gave better results than the classical moments. The classical P ZM and OF M M totally failed to extract the elliptic region because their behavior was too inhomogeneous. Surprisingly, the classical ZM gave very good results. Even though its relative standard deviation was not as good than that of ZM I (see Table 1) the contrast (the mean values difference) between the 2 regions in ZM 2,0 was much higher than in ZM I 2,0 , allowing good regional information. If we compare the results obtained using the four moment-invariant sets we can see that OF M M I showed a relatively better performance in the case of the synthetic data.

We applied our method to the real transrectal prostate image (Fig. 1-b). In this data, the speckle within and outside the prostate was more heterogeneous than that in the synthetic US. Moreover, a urethral catheter was placed to protect the urethra from injury during HIFU ablation, and this catheter projected the classical ultrasonic shadow. However we expected to be able to overcome this problem by interactively setting some well adapted seed points (see because they were less capable of expressing the regional information. This was especially true for P ZM and OF M M that have failed to close the boundary of the prostate. Second, it can be seen that the orthogonal moment-invariant sets provided smoother boundaries than the classical sets, probably because they estimated more homogeneous features. Nevertheless it is difficult to significantly compare the different classes of moments on only the graph cut because the final segmentation result is directly dependent on user interaction (variables seed points and λ value). However, in all the cases the segmentation was easier to tune with the orthogonal moment-invariant sets than with the other sets.

In conclusion, the purpose of this section was only to demonstrate that the region and boundary information obtained by moment-invariant sets on ultrasound image can be directly used in a classical segmentation method. Within the min-cut/max flow graph-cut framework, moments invariants showed their ability to extract texture information from the ultrasound image. Moreover the fact that moment invariants were able to extract boundary information clearly improved the segmentation results.

Conclusion

In this paper, we presented three sets of orthogonal moment invariants that can be used to extract texture features when the texture pattern has a concentric organization or when its scale changes in the image, as with the speckle in ultrasound images, for example. The momentinvariant features extracted on a simulated circular ultrasound image demonstrated that they are able to extract the global distribution of the speckle and that they are little sensitive to speckle orientation or scale changes. These results also showed that some features were sensitive to regional texture information and others enhanced the boundaries between two textured regions. The features sensitive to the boundaries between two textured regions can be seen as texture "gradient" and can thus be used in segmentation techniques where the boundary information is needed, i.e. active contours, deformable models, level sets, etc. In order to demonstrate the usability of moment-invariant features, such as region and boundary information, we introduced this joint information into the min-cut/max-flow graph cut algorithm. The resulting segmentations on simulated and real ultrasound images demonstrated the benefits of the joint region-and boundary-sensitive moment-invariant texture features. Moment invariants are so eligible to be used to adapt techniques based on the combination on regional and boundary information to the segmentation of ultrasound images. "background"). The estimated boundaries using ZM I features with: b) λ = 0, with only the regional information, and c) λ = 1, with a mixture of region and boundary information.

Fig. 8 -

 8 Fig.8-a). The input images were the same moment features (same order and repetition) as for the synthetic data. However we changed the window size to 11 × 11 to adjust it to the transrectal prostate speckle dimension. Fig 8 compare the segmentation results of the mincut/max-flow graph cut algorithm with a λ = 1, using the different moment sets. Because of the non-homogeneity of speckle within the prostate, all the segmentation results show some gaps within the extracted region . However, two facts can be observed from Fig 8. First, the classical moment sets were less efficient than the invariant moment sets for describing the prostate, mainly

5 Acknowledgment

 5 This work was part of the MULTIP project (ANR-09-TCS-011-04)) which was supported by the French National Research Agency under the Health Technology Program (ANR TecSan), by the National Basic Research Program of China under Grant 2011CB707904, the National Natural Science Foundation of China under Grants 61073138, 60911130370, 61271312, and by a PhD grant from the China Scholarship Council.

Figure 2 :

 2 Figure 2: Invariance ability. Variation of P ZM 3,0 , P ZM I 3,0 in relation to: a) speckle orientation and b) speckle scale. In this later case, the speckle remains constant but the moment window size is varying.

Figure 3 :

 3 Figure 3: P ZM I 1,1 values estimated with a 11 × 11, 21 × 21 and 31 × 31 window size. In this figure only a rectangular window around the elliptic region of the simulated image is shown. The elliptic region boundary is shown in red.

Figure 4 :

 4 Figure 4: Influence of the order p. P ZM I p,0 with p increasing from 1 to 3. In this figure only a rectangular window around the elliptic region of the simulated image is shown. The elliptic region boundary is shown in red.

Figure 5 :

 5 Figure 5: Selection of texture features extracted from classical moments and moment invariants at several order p and repetition q. The first column shows Hu 1 and several orthogonal moments and moment invariants of repetition q = 0. The second column shows Hu 2 and several orthogonal moments and moment invariants of repetition q = 1. In this figure only a rectangular window around the elliptic region of the simulated image is shown. The elliptic region boundary is shown in red. 21

Figure 6 :

 6 Figure 6: Some features extracted on the transrectal prostate ultrasound image.

Figure 7 :

 7 Figure 7: Graph cut segmentation results. a) Graph cut seed points (green: 'object", red:

Figure 8 :

 8 Figure 8: Transrectal prostate ultrasound image segmentation. a) the seed points (green: "object", red: "background") and (in black) the contour drawn by an expert. Others: segmentation results.

Table 1 :

 1 Comparisons of the invariance properties of the moments on a simulated US image.

	RSD Relative Standard Deviation.
	Repetition q = 0	Repetition q = 1
	Method	RSD	Method	RSD
	ZM 2,0	28.1%	ZM 1,1	1125%
	ZM I 2,0	0.70%	ZM I 1,1	77.2%
	P ZM 3,0		P ZM 1,1	1125%
		31.5%
	OF F M 3, 0		OF M M 1,1 2604%
	P ZM I 3,0		P ZM I 1,1	382%
		2.36%
	OF F M I3, 0		OF M M I 1,1 137%
	Hu 1	21.3%	Hu 2	96.8%

Table 2 :

 2 Initial weights for the graph cut segmentation.

	links	costs	condition
	{u, v} λ •	ZM I ′ n,1 (u)+ZM I ′ n,1 (v) 2	for {u, v} ∈ N
		∞	for u label. "object"
	{u,S}
	P object (ZM I n,0 (u))	for the other u
		∞	for u label. "backg"
	{u,T}
	P backg (ZM I n,0 (u))	for the other u

Table 3 :

 3 Comparisons of moments-based texture segmentation on a simulated US image

		Fract area diff Sensitivity Accuracy
	Hu	-7.8%	91.81%	91.42%
	ZM	-4.63%	95.0%	94.63%
	ZM I	-6.42%	93.47%	93.36%
	P ZM	-90.87%	9.13%	9.13%
	P ZM I	-4.22%	95.33%	94.88%
	OF M M	-90.83%	9.17%	9.17%
	OF M M I	-2,35%	96,99%	96,33%

S and T are the two terminal nodes associated with the labels "object" and "background" respectively, see[START_REF] Boykov | Graph cuts and efficient N-D image segmentation[END_REF]

List of Figures


Zulma L Sandoval

Jean-Louis Dillenseger

Evaluation of Computed Tomography to Ultrasound 2D Image Registration for Atrial Fibrillation Treatment

This work aims to evaluate intensity-based similarity measures used in rigid or elastic registration techniques to align preoperative computed tomography (CT) and transesophageal ultrasound (US) images of the left atrium and the pulmonary veins, in order to guide high intensity focused ultrasound (HIFU) ablation therapy in patients with atrial fibrillation. The evaluation of the registration techniques has four steps: (i) design of anatomical 2D models including the left atrium and the pulmonary veins from the cryosections of the Visible Human Project, (ii) application of a known rigid/nonrigid transformation in order to produce a gold standard, (iii) creation of a pair of synthetic US and CT images, (iv) estimation of the transformation by rigid/elastic registration of the simulated images and comparison with the gold standard. The performance of eight similarity measures was evaluated for rigid registration: mutual information (MI), normalized mutual information (NMI), entropy correlation coefficient (ECC), joint entropy (H), point similarity measure based on MI (PSMI), energy of histogram (E), correlation ratio (CoR) and Woods criterion (WC). The evaluation protocol estimated the accuracy and the global optimization tendency around the gold standard. For the elastic registration the level of accuracy obtained using MI and NMI was measured.

Introduction

Ablation procedures have proven to be some of the most effective methods in treating atrial fibrillation [START_REF] Cox | Atrial fibrillation II: rationale for surgical treatment[END_REF]. These procedures aim to establish a line of lesions around the pulmonary veins in order to block trigger points of atrial fibrillation. Ultrasound-guided HIFU is a minimally-invasive alternative to other ablation techniques [START_REF] Pichardo | New design for an endoesophageal sector-based array for the treatment of atrial fibrillation: a parametric simulation study[END_REF]. In this procedure, the ablation path is defined in a preoperative stage on 3D cine CT scans. The ablation itself is performed under 2D ultrasound (US) guidance. The 2D-3D registration of intra-operative US images to the pre-operative CT images is then needed in order to transfer and follow the ablation path in the pre-operative context.

The registration procedure of two images of the same scene acquired at different times or from different points of view and/or by different sensors, aims to align them into a common referential. Usually, one image is fixed as a reference and the other is moved and compared with it. The registration procedure usually consists of 1) applying a geometric transformation to the moving image, 2) comparing similarities between the two images and 3) optimizing the geometric transformation parameters in order to maximize the similarity between the two images.

An intensity-based cardiac US to cardiac CT image rigid registration method has been described in [START_REF] Huang | Dynamic 2D ultrasound and 3D CT image registration of the beating heart[END_REF]. The authors of this paper proposed a spatio-temporal registration of the beating heart using mutual information (MI) metric to drive the spatial alignment. This technique has been improved by performing an intensity-based registration on only the relevant features of the US image after median filtering and low intensity thresholding [START_REF] Huang | Rapid registration of multimodal images using a reduced number of voxels[END_REF]. The question of MI validity thus arises. In CT images, tissues are characterized by intensities distribution. However, in US images, each tissue is inherently characterized by a specific spatial distribution of speckle rather than a specific distribution of gray levels. Therefore, the use of an intensity-based SM in a multimodality image registration problem involving US images is still an open question. Beside MI, other intensity-based SMs have also been proposed in the literature. Some of these metrics could be more adequate than MI in driving the registration in all clinical context. Moreover, rigid registration seems not to be the most appropriate in our beating hearth context. The evaluation of SM should also be performed in an elastic registration context. This document is structured as follows: the geometrical transformation and the similarity measures used in the evaluation are described section 2 and section 3, respectively. In section 4, the evaluation protocol is defined. The experimental results are shown and discussed in section 5.

2.

Geometrical transformations

Image registration can be formulated as an optimization problem in which the cost function C is minimized as:

T µ = arg min Tµ C(T µ ; I F , I M ) (1)
where T µ indicates the parametrized geometrical transformation, I F is the fixed image and I M is the moving image.

The vector µ contains the values of the transformation that we want to optimize in order to maximize similarity between the two images. In this paper we will explore two 2D transformations: Rigid (Eq. 2) and Elastic B-splines based (Eq. 3).

T µ = T (φ, d x , d y ) =   cosφ -sinφ d x sinφ cosφ d y 0 0 1   (2)
with φ = rotation, d x = translation in x and d y = translation in y directions.

T µ (x) = x + x k ∈Nx p k β 3 x -x k σ (3)
with x k the control points, β 3 () the cubic multidimensional B-spline polynomial [START_REF] Unser | Splines: A perfect fit for signal and image processing[END_REF], p k the B-spline coefficient vectors, σ the B-spline control points spacing, and N x the set of all control points within the compact support of the B-spline at x.

Similarity measures

The similarity measures are cost functions that measure the similarity between the I F and the transformed I M . Ideally, the SM has an optimum at the set of parameters that best align the moving image, and values that decrease monotonically with the distance to the optimum. Similarity measure is the critical aspect of an intensity based registration. Its variational behavior (e.g. presence of global minima,...) around the global optimum impacts directly the optimization step and the registration quality.

In this paper, eight intensity-based SMs are evaluated. Six of these SMs use the information from the histogram of images while two of them use the spatial information and intensity values. They are: mutual information (MI) [START_REF] Maes | Multimodality image registration by maximization of mutual information[END_REF], normalized mutual information (NMI) [START_REF] Studholme | An overlap invariant entropy measure of 3D medical image alignment[END_REF], entropy correlation coefficient (ECC) [START_REF] Maes | Multimodality image registration by maximization of mutual information[END_REF], joint entropy (H) [START_REF] Maes | Multimodality image registration by maximization of mutual information[END_REF], point similarity measure based on MI (PSMI) [START_REF] Rogelj | Point similarity measures for non-rigid registration of multi-modal data[END_REF], histogram energy (E) [START_REF] Bro-Nielsen | Rigid registration of CT, MR and crytosection images using a GLCM framework[END_REF], correlation ratio (CoR) [START_REF] Roche | The correlation ratio as a new similarity measure for multimodal image registration[END_REF] and Woods criterion (WC) [START_REF] Woods | MRI-PET registration with automated algorithm[END_REF].

Image registration evaluation

The evaluation of the registration techniques has four steps: (i) design of anatomical 2D models including the left atrium and the pulmonary veins from the cryosections of the Visible Human Project (male and female), (ii) application of a known rigid/nonrigid transformation in order to produce a gold standard, (iii) creation of a pair of synthetic US images and CT slices, (iv) estimation of the transformation by rigid/elastic registration of the simulated images and comparison with the gold standard.

Anatomical models and synthetic US and CT images

We create a set of four synthetic phantom images of the left atrium and pulmonary veins (LAPVs). For these, we first pick four slices from the male and female cryosections of the Visible Human Project1 . As an example, Fig. 1-(a) shows a slice from the male dataset. These LAPVs images were segmented and labeled into six types of tissues: air, fat, bone, muscle, water and blood (Fig. 1-(b)).

Synthetic CT images were obtained by applying the Hounsfield value corresponding to the type of tissue modulated by Gaussian noise (Fig. 1-(c)).

For simulating US images we characterized each tissue by its acoustical impedance and a spatial distribution of speckle. This information is set as an input to an US image simulator [START_REF] Dillenseger | Fast simulation of ultrasound images from a CT volume[END_REF] (Fig. 1-(d)).

Rigid registration evaluation framework

Rigid registration SM has been globally evaluated using an optimization independent protocol proposed by Skerl. In this method the golden standard registration transform is known. The main idea of this protocol is to sample some position in a normalized parameter space around the gold standard. The performance of the SM is then estimated by some statistics of the SM computed on these positions. This protocol allows to estimate the accuracy (ACC) of the SM (the distance between the global maximum and the gold standard) and 4 other statistics related to the behavior of the SM during the optimization: the number of local minima (NOM); the capture range (CR) defined as the distance between the global maximum and the closest local minimum; the measure of distinctiveness (DO) defined as the average change of the SM near the global maximum; and the risk of non convergence (RON) which can be seen as the area of positions around the global maximum to which optimization may converge to local maxima.

Elastic registration evaluation framework

Elastic registration is difficult to evaluate because the same moving image can be obtained by applying different deformation fields on the same fixed image. So the direct evaluation by comparing an estimated deformation field to the gold standard one has no sense [START_REF] Klein | Evaluation of optimization methods for nonrigid medical image registration using mutual information and B-splines[END_REF]. We only assess the influence of the SM on the registration accuracy using the following protocol. First, we estimated the deformation field between the simulated US and CT images using our elastic registration method. We applied then the estimated deformation field to the phantom image. Finally, the difference between this deformed phantom model and the gold standard phantom model is measured by computing the overlap of the several structures areas between both images by the Dice score.

Dice = 2|A 1 A 2 | |A 1 | + |A 2 | (4)
with A 1 the structure surface in the gold standard phantom model and A 2 the structure surface in the deformed phantom image.

Results and discussion

Two registration cases has been studied: the registration of the simulated CT images with the corresponding raw US simulated images (denoted RI in the rest of the paper) and the CT images with the US images preprocessed as in [START_REF] Huang | Rapid registration of multimodal images using a reduced number of voxels[END_REF] (denoted PI). The goal of the evaluation is to measure the impact of the several SMs on the rigid and on the elastic registration.

Rigid registration evaluation: The statistical results for the several SM using Skerl's protocol can be seen in Table 1. In this table the best values are bolded.

It can be observed in Table 1 that WC presents a globally better performance than the other SMs. This is especially true for the ACC, CR and DO indicators. MI and CoR equally give globally relatively high performances. Results show that MI can be used in the registration of CT to US images as it has been done in [START_REF] Huang | Rapid registration of multimodal images using a reduced number of voxels[END_REF]. However, results also show that preprocessing the images does not improve the performance. If we compare the behavior between the raw and preprocessed images, it can be concluded that the performance is almost the same for all SMs except for WC which has a optimal ACC on the preprocessed images. All SMs have a small CR, as was expected for these modalities. This result is consistent with that reported in [START_REF] Skerl | A protocol for evaluation of similarity measures for rigid registration[END_REF].

However a global comparison between the several SM is difficult to perform by the single analysis of Table 1. We propose to compute a composite index as a combination of the statistics for each SM. For this we used principal components analysis (PCA) for ranking the similarity measures. The coefficients of the first principal component give the composite index [START_REF] Mishra | On construction of robust composite indices by linear aggregation[END_REF]. The coefficients of Elastic registration evaluation: The elastic registration has been performed using the dedicated software Elastix [START_REF] Klein | Elastix: a toolbox for intensity-based medical image registration[END_REF]. According to our registration problematic this software proposes several transformations including B-splines. But the choice of SMs suited for multimodality registration supported by Elastix is only limited to MI and NMI. With this restriction the average Dice result for all LAPVs images indicates an overlap area of 93.02 percent for MI and of 97.08 percent for NMI using RI US images. With PI US images the results were: 94.17 percent for MI and 97.94 percent for NMI. These high scores confirm the validity of the elastic registration between US and CT images. For elastic registration NMI seems to be more efficient than MI. In addition, preprocessing US images improve the registration. However, the implementation of the others SMs especially WC or CoR into Elastix toolbox should be investigated.

Conclusions

The evaluation of CT to US images registration was performed for rigid and elastic transformations. The performance of different intensity based similarity measures was compared. The results validate the use of intensity measures to perform the registration between CT and US despite intensity representation in the images. For rigid registration, an optimization independent evaluation of eight intensity similarity measures was performed. The similarity measures WC and MI have the global better performance. For elastic registration, NMI show a slightly better accuracy result than MI when the registration was performed in Elastix software.

Figure 1 .

 1 Figure 1. Images used for evaluation. Left atrium and pulmonary veins phantom images: cryosection of male human thorax (a); segmented and labeled into six different types of tissues (b); corresponding simulated CT (c) and US (d) images.

Figure 2 .

 2 Figure 2. Composite index for the several SMs

Table 1 .

 1 Value of the ACC(mm), RON (10 6 /mm), NOM, CR(mm) and DO (1/mm) properties for the eight similarity measures evaluated on: the raw LAPV image (RI) and the preprocessed LAPV image (PI). Best values are bolded.

		ACC		RON	NOM	CR		DO
		RI	PI	RI	PI	RI	PI	RI	PI	RI	PI
	MI	0.17 0.17 5.4	6.4 100 93 0.046 0.046 6.7 8.8
	NMI	0.21 0.19 12.7 13.0 72	58 0.046 0.046 7.5 9.5
	ECC	0.26 0.24 12.2 12.6 89	86 0.057 0.057 7.1 9.3
	H	0.29 0.29 34.2 34.7 38	34 0.046 0.046 4.8 4.3
	PSMI 0.22 0.29 34.1 34.3 50	40 0.046 0.046 1.8 1.6
	E	0.29 0.29 33.2 34.0 46	37 0.046 0.046 5.9 5.2
	CoR	0.14 0.14 6.2	7.4	92	84 0.046 0.057 8.8 11
	WC	0.12	0	7.7 10.4 111 100 0.046 0.079 14 4.4

Computing in Cardiology 2013; 40:245-248.

http : //www.nlm.nih.gov/research/visible/visible gallery.html