Philippe Lasaygues

Serge Mensah

Régine Guillermin

Julien Rouyer

Emilie Franceschini

Non-linear inversion modeling for Ultrasound Computer Tomography: transition from soft to hard tissues imaging

Keywords: Ultrasound Computer Tomography, Inverse Born Approximation, Elliptical Projection, Iterative Approximation, Soft Tissues Imaging, Hard Tissues Imaging, Breast, Bones

Ultrasound Computer Tomography (UCT) is an imaging technique which has proved effective for soft-tissue (breast, liver,…) characterization. More recently, the use of UCT has been envisaged for bone imaging. In this field, the large variations of impedance distribution (high contrast) require that a finer model of wave propagation be integrated into the reconstruction scheme. Here, the tomographic procedure used is adapted to broadband data acquired in scattering configurations while the heterogeneous objects (Born approximation) are probed by spherical waves. An "elliptical" Fourier transform has been derived to solve the near-field inverse problem. This transform differs from the standard Fourier Transform in that, instead of plane waves, families of harmonic ellipsoidal waves are considered. For soft tissues it is possible to separate the impedance and speed of sound contributions and to reconstruct their cartographies using dedicated near-field Radon transforms. In the case of highly heterogeneous media such as bones, iterative inversion schemes are proposed. The various reconstruction procedures are set against experiments.

INTRODUCTION

This paper presents the theoretical framework that enables one, with a unified formulation, to find specific responses, on the one hand, to the problem of soft tissue imaging (dedicated to breast cancer detection) and, on the other hand, to the problem of bone characterization (dedicated to osteoporosis, bone infection and cancer). The difficulties raised are somewhat different as, in soft tissues, the very small fluctuations to be quantified suffer from their very low values. This poor echogenic index generally induces low detection probability, for instance in the case of large diffuse masses. In addition, invasive lesions that must not be overlooked may be of millimetric size. In bone imaging, the difficulties arise from the very high contrast that strongly alters the propagation of the ultrasonic waves. Solutions consist in optimally assessing these non-linear effects in an iterative approach aiming at local linearization. In this article we describe a near field UCT method based on the use of the first-order Born approximation, applied to the case of an homogeneous and constant background. The unknown object function, which is assumed to be weakly heterogeneous, is linearly related to the field measured via a tool called the Elliptical Fourier transform. The inverse problem is related to an extension of the filtered back-projection algorithm. This technique is suitable for breast inspection where the probe is either in contact with the skin or located within a near field distance when using a coupling device (water bag or water tank). However the first-order Born UCT has some limitations when dealing with highly contrasted scatterers. When the problem can be reduced to the study of a fluid-like cavity buried in an elastic cylinder surrounded by water, we have proposed an extension of this method, by taking into account physical phenomena such as wave refraction [START_REF] Ouedraogo | Multi-step compensation technique for ultrasound tomography of bone[END_REF][START_REF] Teng | Application of ultrasonic tomography in industrial materials[END_REF][START_REF] Ouedraogo | Contrast and velocity ultrasonic tomography of long bones[END_REF][START_REF] Lasaygues | Progress toward in vitro quantitative imaging of human femur using Compound Quantitative Ultrasonic Tomography[END_REF]. The approach, named Compound Quantitative Ultrasonic Tomography -CQUT, is purely experimental and consists in performing reflection and transmission measurements, using an iterative correction procedure, which compensates for refraction effect arising at the boundary between bone and the surrounding tissues. The main limitation of the CQUT method is the heavy experimental-costs involved (multiple iterative experiments), and we suggest a purely numerical non-linear inversion algorithm, named the Distorted Born Diffraction Tomography -DBDT, in which the application of the Born approximation is extended to higher orders [START_REF] Lu | Image reconstruction with acoustic measurement using distorted Born iteration method[END_REF]. In comparison with the previous upgrade, DBDT combines diffraction measurements [START_REF] Lasaygues | Ultrasonic Computed Tomography[END_REF]. The performances and limitations of these tomographic methods applied to breast and bone imaging problems are presented.

FORMULATION

Soft-tissue forward problem

Ultrasound Computer Tomography (UCT) applied to soft tissues, (we confine ourselves to breast inspection), has been studied in several scientific research projects [START_REF] Greenleaf | Introduction to Computer Ultrasound Tomography[END_REF][START_REF] Andre | A new consideration of diffraction computed tomography for breast imaging: Studies in phantoms and patients[END_REF] and the use of powerful computers makes it possible nowadays to introduce rather complex algorithms [START_REF] Mensah | Enhanced Compressibility Tomography[END_REF][START_REF] Franceschini | An optimization method for quantitative impedance tomography[END_REF]. Numerous experimental devices have also been developed [START_REF] André | High-speed data acquisition in diffraction tomography system employing large-scale toroidal arrays[END_REF][START_REF] Waag | A ring transducer system for medical ultrasound research[END_REF][START_REF] Rouyer | Conformal Ultrasound Imaging System for Anatomical Breast Inspection[END_REF]. In these approaches, breast tissues are described by means of acoustical models, the mechanical parameter pairs can be either (density, compressibility) or (density, speed of sound); we adopt an (impedance, speed of sound) formulation since these parameters offer separate contributions in transmission and in reflection. Let ρ 0 , c 0 be the acoustic characteristics of the homogeneous coupling medium, and let us introduce the contrast parameter functions of speed of sound and of the acoustical impedance:

α x ( ) = c 2 x ( ) -c 0 2 2 c 2 x ( ) and ξ x ( ) = Log z z 0 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ , (1)
the governing pressure propagation equation is written:

- 1 c 0 2 ∂ 2 p ∂t 2 + ∇ 2 p ≡ S ≡ - 2α c 0 2 ∂ 2 p ∂t 2 -∇α∇p + ∇ξ∇p. (2)
The object to be inspected is excited with a spherical wave having an infinite frequency band:

P i x, e, t ( ) = δ t - x -e c 0 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ 1 4π x -e (3)
Let us consider the differential distribution:

υ e x, t ( ) = -n e 1 x -e δ t ( ) -1 c 0 δ ' t ( ) ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ (4)
unit vector n e is shown on Figure 1, the prime stands for the first derivation and let us note:

S x, e, t ( ) = S * x, e, t ( ) * t p i x, e, t ( ) = - 2α c 0 2 ′′ δ t ( ) -∇α -∇ξ [ ] υ e ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ * t p i x, e, t ( ) , (5)
The temporal support of ( )

t e x S
, , * being of null width (sum of derivatives of Dirac's distribution), in case of spherical transmitted waves (small pitch transducer), one can easily derive the integral expression of the scattered field :

( ) ( ) ⋅ - - π ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - - - - = * - - π ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - - - δ - = ∫ ∫ Σ ℜ t t d x d e x x
Ellipse t Σ of foci ( ) r x, is defined by the relation c 0 t = rx + xe . Then, p d x, e, t ( ) is the elliptical projection onto t Σ of the source terms. Given the two points e and r, the former expression can be rewritten in its integral form marked by « * )

»: G ) * S r, e, t ( ) = δ t -′ t - r -′ x c 0 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ 4π r -x ' S x', e, ′ t ( ) d x' dt ' ℜ 3 ∫ = δ t -′ t - r -′ x c 0 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ 4π r -x ' S * * t ' p i x ', e, t ' ( ) d x'dt ' ℜ 3 ∫ (7)
Hence the synthetic expression of the scattered field (Born approximation) :

p d r, e, t ( ) = G ) * S r, e, t ( ) = G ) * S * * t p i r, e, t ( ) (8)
If instead of a spherical wave of infinite bandwidth one uses a harmonic spherical transmit wave, then [START_REF] Greenleaf | Introduction to Computer Ultrasound Tomography[END_REF] becomes:

p d r, e,ω ( ) = - e ik r-x + x-e ( ) 16π 2 r -x x -e S * x, e,ω ( ) d x Σ t ∫ υe x, e,ω ( ) = n 1 x -e -ik ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ S * x, e,ω ( ) = k 2 2α -∇α -∇ξ [ ] υe x, e,ω ( ) ⎧ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ (9)
The first of these expressions reveals a transform similar in its form to the standard Fourier transform. It is specified by the harmonic decomposition functions having as spatial support a time-indexed ellipsoid, and for foci the loci of the transmitter and of the receiver. This decomposition has been named the Elliptical Fourier transform (EFT). We have shown in a previous paper [START_REF] Mensah | Near field tomography[END_REF] that this transform is invertible. Thus, the scattered field is the ellipsoidal spectrum of the secondary source term specified by the acoustical parameters. Nota 1: In the harmonic regime, since the harmonic functions are the eigen functions of the temporal convolution operator, (8) is written:

p d r, e,ω ( ) = G ) * S r, e,ω ( ) = G ) * S * × p i ( ) ⎡ ⎣ ⎤ ⎦ r, e,ω ( ) (10)
Nota 2: When no confusion is possible, for any well-behaved function f, ( ) ω , , e r f refers to the spectrum of function f r, e, t ( ) . When distinction is required the spectrum is also noted: f r, e,ω ( ) = f r, e,ω ( ) = f = F f r, e, t ( )

⎡ ⎣ ⎤ ⎦ , F is
the standard Fourier transform. Nota 3: When source term S * corresponds to a harmonic plane wave, excitation can be easily detailed from [START_REF] Mensah | Enhanced Compressibility Tomography[END_REF]. Indeed, in that case, the source locus e is "pushed away" to infinity and υe x, e,ω ( ) = -n 0 ik (n 0 is the unit vector in the incident direction). Thus, [START_REF] André | High-speed data acquisition in diffraction tomography system employing large-scale toroidal arrays[END_REF] This is the result obtained in ultrasound tomography development using far-field measurements and standard Back-Projection (BP) algorithm [START_REF] André | High-speed data acquisition in diffraction tomography system employing large-scale toroidal arrays[END_REF]. We observe here the possibility of continuously adjusting the inversion operator from far field down to near field thanks to υe x, e,ω ( ) vector.

( ) ( ) ( ) ( ) [ ] ω ξ ∇ - ω α ∇ - ω α = ω , , , 2 , 0 2
* x x n ik x k x S
Then, applying the EFT (overlined with an arrow) to the formulations in [START_REF] Andre | A new consideration of diffraction computed tomography for breast imaging: Studies in phantoms and patients[END_REF] and considering the specific derivation rules related to the elliptical analysis [START_REF] Mensah | Near field tomography[END_REF], one obtains:

υ x,ω ( ) = υe x,ω ( ) -υr x,ω ( ) K = kn ϕ,θ -p d K ( ) = k 2 2α sr K ( ) -υe ⋅ υ s r uuuuu ( ) * α sr -ξ sr ⎡ ⎣ ⎤ ⎦ K ( ) ⎧ ⎨ ⎪ ⎪ ⎩ ⎪ ⎪ (12)
The scattered field is the result of the convolution in the elliptical spectral domain of the elliptical spectra with the directivity index functions. If the transmitter (resp. receiver) /pixel distances are of the order of the wavelength, we obtain with a good approximation

υe ⋅ υ ≈ -k 2 1+ n e .n r ( ) , where n l = ∇ x -l ( ) (13)
e n is the unit vector in the direction transmitter (resp. receiver) / voxel. Then, the system is reduced to:

p d K ( ) = -k 2 1-n e . n r s r uuuuuuu ( ) * α sr -1+ n e . n r s r uuuuuuu ( ) * ξ t ⎡ ⎣ ⎤ ⎦ K ( ) (14)
Indeed, these equations involve two acoustical parameters offering distinct directivity patterns that can be separated according to the spatial scanning performed:

• n e = -n r , p d = k 2 α sr , in pure transmission, the scattered field is conditioned by the quadratic fluctuations of the sound speed.

• n e = n r , p d = k 2 ξ t , in pure reflection, the scattered field results from impedance fluctuations. These results are coherent with those obtained using plane wave excitation and plane wave decomposition (standard Fourier transform). Here, they are extended in the sense that the analysis has been done taking into account the curvature of the "elliptical" wavefront.

"Near-field" tomography

By applying the inverse elliptical Fourier transform to the elliptical spectra, it is possible to reconstruct the impedance or the sound speed maps since they are respectively obtained from pure reflection and transmission measurements.

Fast inversion procedure

For this purpose, the acquisition procedure consists in filling in the respective elliptical Fourier domains with the highest sampling densities by performing both a spatial scan and a frequency sweep using broadband signals. In practice, however, the transducers have a finite bandwidth and one may be able to reconstruct only a bandpass (respectively, low-pass) filtered version of the impedance (respectively sound speed). Besides, rather than using multidimensional inverse elliptical Fourier transforms to retrieve the object from diffraction measures, it is better to use a reconstruction procedure based on classical one-dimensional (1D) Fourier transforms that are numerically optimised. Therefore, we have defined a near-field extension of the Radon transform and introduced an elliptical backprojection algorithm similar in its form to the standard backprojection (FBP) algorithm.

Let f x ( ) be a well-behaved object function defined over the compact domain D of 2 . f represents one of the separable parameters specified above. For a given configuration angle θ = θ 0 , (cf. Figure 1), we define the elliptical Radon transform for any pair (s, φ) ∈ {s ∈ , 0≤φ≤π} by:

R ε f [ ] s,ϕ ( ) = f x ( ) δ s -r -x -x -e ( ) 16π 2 r -x . x -e ∫ d x , (15)
and consider the projection p φ (s) under incidence φ:

R ε f [ ] s,ϕ ( ) ϕ fixed = p ϕ s ( ) (16)
We introduce the scalar product:

n ϕ,θ •x = -r -x + x -e ( ) . (17)
Here, its negative value results from the fact that the wave is propagating toward the origin, and is opposed to the outgoing orientation of axis x.

We derive the synthesis of function f from its ellipsoidal decomposition (projections p φ (s) under incidence φ): [ ] ϕ = 16π 2 r -x x -e , accounting for the divergence of the beam, one can retrieve all the constitutive elements.

f x ( ) = 1 2π ( ) 2 pϕ S ( ) -∞ ∞ ∫ 0 π ∫ S e iS ϕ •x x [ ] ϕ dS dϕ
A difference however, which is worth noticing: the nature of the projection differs since, here, we consider an elliptical projection p φ (s) over Σ t defined by:

s = c 0 t = S ϕ •x = r -x + x -e , (19)

Elliptical Backprojection (EBP) operator

Let us consider an arbitrary function h(s,φ) where s = n ϕ • x , the elliptical backprojection operator is:

B h s,ϕ ( ) ⎡ ⎣ ⎤ ⎦ x ( ) = 1 2 × 1 π x [ ] ϕ 0 π ∫ h n ϕ •x,ϕ ( ) dϕ. ( 20
)
Function f is then recovered by EBP [START_REF] Waag | A ring transducer system for medical ultrasound research[END_REF]:

f = B o F -1 o Abs o F o R ε o f = R ε -1 o R ε o f (21)
where the response of filter Abs is Abs(S) = |S|. Let us call p dm the "measured scattered field" and S est * the estimated parameter map reconstructed from the data acquired; remember that S est * is bound either to the impedance map or to the speed of sound according to the protocol of acquisition (pure reflection or transmission). Under such conditions,

S est * = R ε -1 p dm [ ] (22)
Thus, the Elliptical Backpropagation Algorithm leads to an Inverse Born Approximation (IBA) with a (almost) "constant background" (considering the matching of the mean acoustical impedances).

"Hard-tissue" tomography Modified elliptical Radon transform:

With hard biological tissues having larger acoustical impedances and speed of sound than those of the surrounding medium, the weak scattering assumption is no longer realistic. Indeed, the transmitted wave is "deviated" from its initial course because of refraction. This results in a new distribution of the energy initially transmitted to the Region Of Interest (ROI). Consequently, only a small percentage of the incident energy penetrates into the contrasted medium, and the recorded signal is highly attenuated, especially as far as the high frequencies are concerned. Thus, if it is not taken into consideration prior to the inversion, this redistribution may lead to a loss of information; for instance, strong geometrical distortions or large shadows may occur. On the contrary, appropriate scanning of the highly contrasted medium can provide good qualitative imaging. In addition, if prior information is available, quantitative recovery of the parameter distributions may be expected. Thus, due to the important sound speed and impedance fluctuations, the assumption that sound speed is constant in relations ( 7) and ( 8) is no longer valid. Instead, one has to take into account the spatial dependence of the speed of sound in the derivation of the scattered field:

p d r, e, t ( ) = - S * x, e, t + τ x ( ) - r -x c 0 - x -e c 0 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ 16π 2 r -x x -e d x Σ t ∫ (23)
Here, the temporal correction is approximated on the basis of straight line propagation inside the hard medium:

τ r, e, x ( ) = 1 c 0 c x' ( ) -c 0 c 0 d x' ex + xr ∫ = τ e x ( ) + τ r x ( ) (24)
Hence, the overall correction is the sum of the respective temporal corrections along the emitter-voxel (located in x) path and along the voxel-receiver path. Indeed, this correction term is unknown since the sound speed parameter is directly or indirectly sought for; it is generally estimated via an iterative tomographic procedure initialized on the basis of prior information that, for instance, may correspond to the bone mean speed of sound for equivalent age and health status of the patient.

Then, the expression of the corrected Radon transform is:

R ε c f ⎡ ⎣ ⎤ ⎦ s,ϕ ( ) = f x ( ) δ s + c 0 τ x ( ) -r -x -x -e ( ) 16π 2 r -x . x -e ∫ d x (25)
The corrected inverse elliptical Radon transform is built from the corrected elliptical backprojection operator:

R ε c ⎡ ⎣ ⎤ ⎦ -1 = B c o F -1 o Abs o F where B c h x ( ) = B h n ϕ • x + c 0 τ x ( ) ,ϕ ( ) ⎡ ⎣ ⎤ ⎦ (26)
The imaging procedure: Contrast imaging, for instance, assessing "cortical thickness" in bone imaging, consists mainly in assimilating the background comprising two areas: the reference surround medium (generally water) and the homogeneous solid part (cortical bone). The default sought for corresponds to a perturbation in that cortical area that creates a cavity in it.

A first practical approach in operating with reflection measurements (reflectivity imaging) consists in using low frequency (≤ 1 MHz) activation in order to take advantage of the large penetration depth. In that case, the Born Approximation (BA) is still met since firstly the wavelength is large compared with the size of the heterogeneity and secondly, the fluctuations of the mechanical parameters inside the cavity are small. The "constant background" inversion procedure can be based on an Elliptical Back-Projection algorithm (20) complemented with specific signal processing refinements (deconvolution). There, despite the artifacts and biases affecting the assessment of the shell thickness, good quality (contrast and resolution) images may be obtained [ 15 ]. In a second approach, one can adopt the Born Approximation (BA) with "a variable background" which consists in considering the hard medium as homogeneous while its geometry and the mean speed of sound are unknown. This requires an iterative inversion methodology, which is generally used whenever quantitative information (speed of sound distribution) is sought for. Indeed, in this case, finding solutions involves either using non-linear schemes or performing extensive studies that expand the boundaries of the use of the Born series. The so-called Distorted Born Iterative (DBI) method is a mixed approach that implicitly specifies both the field inside the object and the Green function. Then starting from [START_REF] Andre | A new consideration of diffraction computed tomography for breast imaging: Studies in phantoms and patients[END_REF], one can write the expression of the measured scattered field:

p dm r,e,t ( ) = G b ) * S * * t p b r,e,t ( ) (27)
G b is the Green function of the "variable background" (the hard medium is considered as homogeneous, but its geometry and its mean speed of sound are unknown), and p b is the associated field so that the total measured field is:

p tm r, e, t ( ) = p i r, e, t ( ) + G b ) * S * * t p b r, e, t ( ) . (28)
The three unknowns in (28) are the Green function, the field inside the object and contrast function S * . These variables are assessed using the following iterative algorithm:

S est *n = R ε c n ⎡ ⎣ ⎤ ⎦ -1 p dm ( ) S est *n x ( ) = c n x ( ) , ξ n x ( ) ( ) τ n-1 r,e, x ( ) = 1 c 0 c n-1 x ( ) -c 0 c 0 dx, τ 0 x ( ) = 0 ex+ xr ∫ R ε c n ⎡ ⎣ ⎤ ⎦ -1 = B c n-1 o F -1 o Abs o F ⎧ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ p b n-1 r,e,t ( ) = p i r,e,t +τ e n-1 r ( ) ( ) b n-1 r,e,t ( ) G r,e,t +τ e n-1 r ( ) ( ) (29)
Nota 1: When the measurements are made in transmission, S est *n x ( ) gives access to the sound speed distribution, c n x ( ) . In reflection, a band pass filtered version of the log-impedance, ξ n x ( ) , may be obtained.

Nota 2:

The last two relations are implicitly obtained but not evaluated.

To begin with, we proposed the Compound Quantitative Ultrasonic Tomography (CQUT), where refraction is taken into account in order to impose (almost) straight ray propagation inside the shell. From a priori knowledge of the geometrical properties and of the average acoustical parameters (speed of sound), the incident and refracted angles are estimated using the Snell-Descartes laws. Then, a compensation procedure is introduced into the scanning (e.g., by using a beam former) in order to obtain parallel paths within the "homogeneous" contrasted medium [START_REF] Lasaygues | Progress toward in vitro quantitative imaging of human femur using Compound Quantitative Ultrasonic Tomography[END_REF]. In an iterative experimental sequence, reflection and transmission measurements enable one to delineate the boundaries of the shell (the cortical region) and to assess the values of the speed of sound along the wave paths. The estimated speed of sound map obtained at iteration n is used as the geometrical and parametrical a priori information for step n+1. The stop criterion corresponds to a difference between the mean velocities calculated for two successive iterations lower than 5 m/s. A second non-linear inversion method (i.e. the Distorted Born Diffraction Tomography (DBDT), which is also based on the DBI method) was therefore investigated. The medium is modeled without any a priori knowledge by performing a simple geometrical discretization of the object. The algorithm involves successive linearizations of the integral representation. The initial guess in the iterative process is provided by the first-order Born approximation. If the solution is known with the order (n-1), the n-order solution will satisfy [START_REF] Haddadin | Ultrasound inverse scattering for tomographic imaging self-focusing array[END_REF]. At each iteration, the algorithm solves a forward diffraction problem in order to calculate the appropriate inhomogeneous Green function G b n-1 and the internal pressure

field p b n-1 .
Contrary to what occurs with the CQUT, the DBDT requires only a single series of experimental diffraction data. However, it is computational time consuming because it involves inversion of huge, full and complex matrix. The matrix inversion procedure is the key point with this method. Generally non-square ill-conditioned matrixes have to be inverted. A mean-square solution can be calculated using a conjugated-gradient method associated with a regularization procedure [START_REF] Lasaygues | Distorted Born diffraction tomography applied to inverting ultrasonic field scattered by noncircular infinite elastic tube[END_REF][START_REF] Guillermin | Quantitative non-linear ultrasonic imaging highimpedance contrast targets -An experimental investigation[END_REF]. To make use of the broadband frequency content of the impulse signal used, the idea is to begin at the low frequencies, which carry overall information, and to gradually inject the high frequencies to simultaneously improve both the qualitative aspects (the resolution) and the quantitative aspects (the parameterization).

EXPERIMENTAL IMAGING

The feasibility of UCT methods was tested on data obtained using ultrasonic scanners (Figure 2), which permits diffraction, reflection and/or transmission measurements (except for DBDT-see §1.1.6). The UCT reconstruction algorithm and the signal processing modules are implemented on a standard personal computer. 1.1.4. Soft-tissue experimental imaging (pure fluid approximation). The breast training phantom CIRS Model 052A mimics the ultrasonic and anatomic characteristics of tissues found on an average patient in a supine position. The flesh-like phantom made of Zerdine® contains cystic and solid masses which appear respectively dark and bright on the ultrasound and tomographic images. In Figure 3, the echographic view shows four inclusions while the tomographic reconstruction (Figure 4) reveals eight inclusions (five cysts, three masses). These differences are due, on the one hand, to a sharper directivity (thinner cross-section) of the conventional probe, which operates at a central frequency more than twice that used by the tomographic array. One the other hand, the scanner performs an entire 2D circular scan (Ø = 30 cm, 720 reflection measurements) around the breast phantom whereas a unique position of the probe (one incident direction) is considered with the echographic system. Comments: With sonographic semiology, cysts are generally the easiest ones to diagnose reliably since they are round or oval anechoic lesions with smooth wall contours and dorsal echo enhancement. On the contrary, on echographic view, the dense solid mass generates bright echoes at the sound-entrance border and a shadow behind the lesion. These two artifacts (echo enhancement and shadowing) are not reproduced with the tomographic approach (Figure 4). Besides, the Zerdine matrix is depicted with important speckle noise on the echographic image. This noise is strongly reduced on the tomographic representation since the point-spread function of the scanner is recognized as being more isotropic and spatially concentrated than that of the echographic device. Furthermore, circular rings are visible round certain cysts (e.g. the larger one); the corresponding echoes are observed on the sinogram (Figure 4 -left). These rings are slightly perceptible in the conventional ultrasound image depicted in Figure 5 -left, where a zoom is performed on the largest cyst (equatorial plane). These views confirm that, within the matrix, an interface exists which delineates a Zerdine layer around the cystic inclusions. Considering the tomographic approach (Figure 5 -right), this "weak" interface is revealed distinctly by the coherent summation performed over a large aperture of the interrogating wavefront. On the contrary, using conventional ultrasounds, the continuity of the (Zerdine) material induces low level backscattering of the focused transmit waves. However, one reservation must be made concerning tomographic representations: the contour of the dense mass appears slightly speculated. This is because of its strong echogenicity, which associates high values with the elliptical backprojections crossing the lesion domain. This effect may be reduced by using a large aperture phased array in reflection that attenuates the reconstruction noise, as in [START_REF] Mensah | Enhanced Compressibility Tomography[END_REF][START_REF] Rouyer | Conformal Ultrasound Imaging System for Anatomical Breast Inspection[END_REF].

1.1.5. Bone experimental qualitative imaging (fluid-like approximation). 3-D UCTs (Figure 6) of child's fibula were obtained from diffracted measurements. The nominal frequency of the transducers was 1 MHz. The resolution and the size of the image were 0.75 x 0.75 mm and 512 x 512 pixels, respectively. The diffracted sinograms included 2048 projections, each including 8192 samples, involving 32 angular transmitter positions combined with 64 angular receiver positions covering an angle of 360°. The sampling frequency used was 40 MHz. A signal-processing tool [START_REF] Lasaygues | Assessing the cortical thickness of long bone shafts in children, using two-dimensional ultrasonic diffraction tomography[END_REF] was used to determine cortical thickness. This tool was based on a segmentation of the final image and a size correction of the inner boundaries, using a priori information on physical parameters (in this case, mean compressional wave velocity in bone = 3500 ± 100 m.s -1 , bone mass density = 1700 kg.m - 3 ). The 3-D UCT images were obtained by superimposing sequential 2-D UCT images (20 cross-sections). The step between two cross-sections was 1 mm. The interpolation scheme used for a given 3-D image was the shape-based method (Matlab®, MathWorks™). Figure 7 shows 2-D UCT of artificial paired bones (fibula and femur). The both bones show up particularly clearly. But, the zone between the thighbone and the fibula is not clearly imaged and it is impossible to measure the distance between the both bones. The copper rod, placed in-between, isn't visible too. This was mainly due to the multiple diffusion, refraction and others phenomenon, occurring. However, these results obtained with human specimens show that it is possible to characterize the size and shape of bones using Ultrasonic Computer Tomography. Focusing on well-contrasted images obtained using UCT show that this method provides efficient means of assessing the cortical thickness, which is known to be an important indicator of bone strength.

"Bone" experimental quantitative imaging (fluid-like approximation).

The aim here is to obtain information about the shape, dimensions and sound speed profile of the studied object. Classical diffraction tomographic methods based on the first-order Born approximation give good results for weakly scattering targets but fail when it comes to give a quantitative information, such as the sound speed profile, for high impedance contrast targets such as bones. In this case, the first-order Born approximation is no longer valid. We propose to use the DBDT to overcome this difficulty. This is an iterative process that successively solves a linearized inverse problem involving the minimization of the discrepancy between the measured and computed scattered fields [START_REF] Lasaygues | Ultrasonic Computed Tomography[END_REF][START_REF] Lasaygues | Distorted Born diffraction tomography applied to inverting ultrasonic field scattered by noncircular infinite elastic tube[END_REF]. The stability of the inversion process is enhanced by means of Tikhonov regularization. An experimental study was performed with a non-circular homogeneous isotropic shell made of artificial resin (NEUKADUR ProtoCast 113 TM ). The mass density of the resin was 1150 kg.m -3 , and the mean velocity of the compressional wave was 2400 + 50 m/s. The surrounding fluid-like medium (and the hollow area) was water at a temperature of 18 ° (mass density: 1000 kg.m -3 , sound speed: 1480 m/s). The transmitter and the receiver were placed at 17.5 cm to the center of the bench. Ultrasounds were generated using a pulser transmitter/receiver device (Panametrics HP 5055PR) and wide-band Panametrics transducers. We used successively two pairs of transducers with a nominal frequency of 0.250 MHz and 1 MHz, respectively. The usable bandwidth range from approximately 135 kHz to 375 kHz and from 500 kHz to 1.5 MHz. The object was placed in the center of the measuring system. The sector scanned was 360 degrees around the target with both transducers with an angular increment of 10 degrees (36 x 36 signals for each paired transducer). Iterations were performed by gradually increasing the working frequency with 3 frequencies chosen in the transducers bandwidth: 150 kHz, 350 kHz and 1 MHz. For iterations at 150 kHz, the fist-order Born approximation solution is used as an initial guess. At other frequencies, the initial guess was the final result of the iteration process performed at the previous frequency (for example, the result of the inversions performed at 150 kHz was used as the initial guess for the inversions performed at 350 kHz, etc). Figure 8 show the first-order Born solution and three sequential iterations of the DBDT. It can be seen that the resolution and the quality of the contrast is gradually improved. The final result of the iteration process was satisfactory. The geometry was fairly accurate (position and dimensions reconstructed with mean relative error of the order of 5%) and the velocity was estimated with a relative error of about 10 %.

CONCLUSION

Ultrasound Computer Tomography (UCT) appears as an alternative imaging tool able to reveal the internal structure of soft tissues (mammography), to delineate the shape of bones (cortical thickness) and even to provide the possibility of parametric estimation (the sound speed map). The theoretical framework corresponds to the analytical formulation (Green function) of the forward problem considering point-like transmitters and receivers. Thus the shape of the transmitted wave (the incident field) is spherical, as are also assumed to be the secondary sources induced by the elementary scatterers of the investigated media.

As regards soft tissue imaging, the derivation of the near-field forward problem has led us to define new elliptical Fourier transforms and associated elliptical Radon transforms. Using these two transform pairs (direct and inverse transforms) the contributions in the scattered field of impedance and of sound speed respectively can be separated. This opens the way to single parametric reconstruction (imaging) i.e. to assessing the distribution of the acoustical parameters using optimized numerical procedures (Fast Fourier Transforms).

For highly heterogeneous media having greater acoustical impedances and speed of sound than those of the surrounding medium, additional consideration of a "variable (homogeneous) background" is necessary if one seeks parametric imaging. Here, the strong non-linearities in the parametric distributions have been dealt with using an iterative inversion procedure.

The examples dealt with in this paper were industrial breast phantoms and real thighbones. In comparison with conventional ultrasounds, UCT methods applied to soft tissues offer reduced speckle noise, artifact-free representation, and enhanced quality (point resolution and contrast) imagery. However, in the vicinity of hard lesions, the presence of strong scatterers may induce some reconstruction noise that could affect the image interpretation (spiculated artifact).

Figure 1 :

 1 Figure 1: Scattering acquisition configuration in acoustic tomography: the punctual source located in e corresponds to the focus of the spherical transmitter, the receiver focus is in r.

( 18 )

 18 This relation resembles the continuous expression of the classical backprojection algorithm. Except the normalization term, x

Figure 2 :

 2 Figure 2: (Left) Soft tissues imaging: a semi-circular multi-channel 3 MHz-antenna, using 1024 piezo-composite elements. The frontend electronics incorporate 32 independent parallel transmit/receive channels and a 32-to-1024 multiplexer unit for fast acquisition of signals (7 sec.) [13] (Right) Hard tissues imaging: a 2D-ring antenna using 8 piezo-composite 1 MHz-transducers, and a 3 MHz-array using 128 piezo-composite elements The front-end electronics incorporate 8 independent parallel transmit/receive channels and a 8to-128 multiplexer unit for fast acquisition of signals (1 sec.)

Figure 3 :

 3 Figure 3: Echographic view (Mindray, 7.5 MHz) of breast training phantom CIRS Model 052A

Figure 5 :

 5 Figure 5: Comparison of conventional (left) and UCT (right) image of the cyst phantom (zoom of the Figure 4 -right).

Figure 6 :

 6 Figure 6: Cortical thickness of a 12-year old child's fibula (a) sample embedded in an agar-gel shell (b) Twenty sequential crosssections were performed in this case with a 1-mm step (c) 3-D tomogram (300 x 300 pixels) (d) the mean dimensions (after processing) of the bone sample were 7 ± 2 mm on the outside and 3 ± 2 mm in the inner cavity.

Figure 7 :

 7 Figure 7: 2-D UCT of adjacent fibula and femur bones, 300 x 300 pixels

Figure 8 :

 8 Figure 8: (a) NEUKADUR ProtoCast 113TM resin target. Results of DBDT: (b) fist-order Born solution at 150 kHz, (c) after iterations at 150 kHz, (d) after iterations at 350 kHz and (e) after iterations at 1 MHz.

Proc. of SPIE Vol. 8320 832004-1

Proc. of SPIE Vol. 8320 832004-2

Proc. of SPIE Vol. 8320 832004-3

Proc. of SPIE Vol. 8320 832004-4

Proc. of SPIE Vol. 8320 832004-5

Proc. of SPIE Vol. 8320 832004-6

Proc. of SPIE Vol. 8320 832004-7

Proc. of SPIE Vol. 8320 832004-8

Proc. of SPIE Vol. 8320 832004-9

Proc. of SPIE Vol. 8320 832004-10

Proc. of SPIE Vol. 8320 832004-11

Proc. of SPIE Vol. 8320 832004-12

Proc. of SPIE Vol. 8320 832004-13

Proc. of SPIE Vol. 8320 832004-14

ACKNOWLEDGMENT

The authors are grateful for medical assistance to Drs D. Amy, radiologist, Aix en Provence, Dr P. Petit and Dr J.-L. Jouve at Public Assistance Hospitals and the "Timone" Children's Hospital in Marseille. Our thanks also go to Eric Debieu, Stephan Devic, and Alain Busso for their technical assistance. This study was based on research supported by the National Center for Scientific Research (CNRS), the Canceropôle PACA, the Provence-Alpes-Côtes d'Azur Council, and the French National Research Agency (ANR).

Theoretical considerations have demonstrated that large aperture measurements in reflection improve the robustness of the reconstruction. This approach is currently under development [START_REF] Rouyer | Conformal Ultrasound Imaging System for Anatomical Breast Inspection[END_REF]. Applied to bone imaging, UCT methods give geometrical information of poorer quality (coarser resolution, geometrical distortions), but the reconstruction obtained is still clinically relevant, especially in the case of children's bones. Besides, the mechanical assessments provided are in direct connection with medical demand (i.e. modeling and strength status of the bone) and it is a non-ionizing modality that can be used in pediatric care. However, estimations of geometry and mechanical information must be associated in a common procedure. Ultimately, the improvement in image quality and in accuracy will be obtained only with a more refined characterization of the 3-D cortical structure. The definition of this structure will specify the "distorted" 3-D Green function, i.e. the conditions of propagation in that cortical structure.


Andrés Arciniegas

Loïc Brancheriau

Philippe Gallet

Philippe Lasaygues

Numerical modeling in quantitative ultrasonic tomography of standing trees

Introduction

The upward-growth structures on the standing trees are changing fast in order to compete for environmental strength and the natural weather forces, and obviously relatively to the age. These changes result in complex mechanical variations, and the physical properties result of gradual transition from juvenile to mature stage into the tree. In addition to naturally occurring material property variations, abnormal diseases content also varies from the pith to the bark. Anyway, healthy trees can present mechanical weaknesses, and therefore be dangerous in urban environments. Studies have demonstrated the effectiveness of acoustics in detecting/imaging trees' decay. The main commercial acoustic devices for tree assessment are essentially based on methods using low frequency acoustic waves. Three tools are "PICUS R Sonic Tomograph" (Germany) [START_REF] Rust | A new tomographic device for the nondestructive testing of trees[END_REF], the "ARBOTOM R " [START_REF] Rinn | Holzanatomische grundlagen der schalltomographie an baumen[END_REF] (Germany), and the "FAKOPP" [START_REF] Divos | Resolution of stress wave based acoustic tomography[END_REF] (Hungary). They use an impulse hammer to excite the trunk and measure the acceleration on circumferential points (accelerometer) with a limited number of sensors. These methods give average values recorded along the acoustic pathway, but they are not always suitable to access 2D resolution (about few centimeters for defect size).

During the last decade, ultrasonic tomography made considerable progress in medical and industrial domain. With the evolution of the electro-acoustical technology new devices for image acquisition have appeared. The main aim of our project is to develop a mechanical ultrasonic imaging device for standing trees (ARB'UST) using two transducers (from 50 to 90 kHz), that could be transportable and used in an in situ condition. With ARB'UST, data is collected sending energy through the tree as many different angles as used to compute an image of the changes of a physical quantity in a slice. The ultrasonic waves cross from the emitter to the receiver at appropriate locations around the tree. These locations are chosen so that as many rays pass through as much of the object as possible. The parameter measured is the time-of-flight (travel time) of the ultrasonic wave between the emitter and the receiver. The final ultrasonic tomogram represents the slowness (the velocity) distribution (related to the axial and lateral resolution of the transducers) in a thin cross-section of the tree (related to the elevational resolution of the transducers).

In ultrasonic tomography different approaches have already been proposed in literature [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF][START_REF] Grangeat | La tomographie : fondements mathématiques, imagerie microscopique et imagerie industrielle[END_REF][START_REF] Lefebvre | Materials and Acoustics Handbook, ch. Acoustic Tomography, Ultrasonic Tomography[END_REF]. There are two major classes of reconstruction techniques. One class is based on the projection-slice theorem such as the filtered backprojection method (FBP). This method is fast and the data must be acquired on spaced sets of straight ray paths (projections). The main actual challenge in tree imaging is the choice of the effective method to inverse the scattered data, and the first step in this work is to test two tomographic inversion methods. This paper focus in an approach by numerical simulation of the reconstructions of cylindrical homogeneous isotropic medium with defects.

Materials and methods

In this paper a comprehensive study is conducted on the reconstruction capability of travel time tomography method for a ring transducer scanner using computed-simulated data for a numerical tree phantom. The scheme to simulate ultrasonicwave propagation through the phantom depends of the experimental set-up and assumptions.

From the mechanical point of view, ARB'UST device was designed to perform sectorial scanning such as fan beam or transmission mode. The general architecture of the system is that of an ultrasound tomographic scanner. An aluminum ring holds two 50 kHz -transducers allowing the bi-static 2D-survey around the trunk [START_REF] Brancheriau | Ultrasonic imaging defects in standing trees -development of an automatic device for plantations[END_REF].

The numerical simulation can be used to study the accuracy of different travel time tomography for different tree phantoms (shapes and physical properties mimics). Numerical tree phantoms are placed inside a ring where transducers are positioned at the boundaries at regular angular spacing locations, from 0 o to 360 o . A 120 o circular sector is not scanned due to positions where the receiver cannot be experimentally placed (Figure 1). A typical angular step of 4 o leads to 61 points per projection. For a complete scan, the total number of acquisitions is therefore 5490. For one projection, the emitter remains in its position while the receiver moves around the ring (Figure 1).

Figure 1: The principle of acquisition

For this purpose, pure compression ultrasonic wave were transmitted in the sample, and no shear waves were taken into account. The medium was supposed homogeneous and isotropic. The wave velocities were independent of frequency. The ultrasonic wave attenuation would be neglected. Only the propagation path phenomena were considered, and the wave delays studied. The numerical tree phantom was then represented by the ultrasonic slowness distribution, and each pixel had a constant value of slowness of the ultrasonic wave. The implementation of the phantom, algorithms and simulations was made in MATLAB R .

Mimicking phantom of tree

Simulations can be easily performed for a large number of well-controlled mimicking phantoms of tree, and can help to understand the features of wavefields acquired by transducers. The simulation can then be used to study the accuracy of travel time tomography for different sizes and shapes of phantom. In this paper, it is proposed to apply this method for a phantom having two main perturbation sizes and shapes inside the host medium and including a small defect. The numerical phantom is a cylindrical homogeneous isotropic object, with cylindrical homogeneous isotropic pertubations (Figure 2). Each small point at the boundaries of the object represent an angular position where the emitter is going to be placed. Table 1 shows the parameters used to construct the phantom.

Forward problem

To test the reconstruction algorithms, it was carried out a numerical implementation of the experimental protocol (acquisition sequences). Equation 1 is known as the Radon transform of a function f(x, y).

P θ (τ) = (θ,τ) line f (x, y)ds (1)
In the reality, the function f represents the physical parameter to be measured. For a fan beam geometry, an independent measurement, or point of projection, is defined as the ray integral (along the line τ) as is given by P θ (τ) for a constant observation angle θ [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]. A projection is therefore the collection of the ray integrals measured along the fan.

Travel time t between an emitter n and a receiver m is represented by a line integral of the slowness p along the ray path Γ [START_REF] Socco | Feasibility of ultrasonic tomography for non-destructive testing of decay on living trees[END_REF], as shown in Eq. 2.

t(n, m) = Γ(n,m) p(x, y)dΓ (2)
The discrete forward problem for an independent travel time measurement is then approximated by a weighted summation for N pixels crossed by a given ray path (Eq. 3).

t(n, m) = N q=1 Γ(n,m) p q d q (3)
The distance D for a given path is assumed as a multiple of constant discrete step d q (pixel), thus D = Nd q . Dividing each side of Eq. 3 by D, the numerical implementation of Radon transform for the slowness P(n,m) is then:

P(n, m) = 1 N Γ(n,m) p(x, y) (4)

Reconstruction by the filtered backprojection method (reference method)

Let us consider a mesh of pixels. To find the value of the slowness at any point of the object, it is assigned the value of the projection at any pixel placed in the path that it covers, and then adding all contributions from all projections. This principle is called backprojection. This method, however, leads to biased values and it is necessary to correct each projection by a specific filtering. The filtered backprojection consists therefore initially to filter each projection using a ramp filter (often combined with a low pass filter to avoid noise amplification), and then to project back the filtered projections for the different observation angles [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF].

Reconstruction by the Layer Stripping based method

In this paper, it is proposed an application of a new inversion algorithm based on a "Layer Stripping" method [START_REF] Franceschini | Tomographie ultrasonore dédiée à la détection du cancer du sein[END_REF]. This procedure is an iterative approach, for transmission measurements in fan-beam mode, where ray paths are analyzed as if the object was discretized by concentric layers. The tomographic reconstruction consists of the local estimation of slowness of an identified concentric layer and to correct it with respect to the values of slowness that were already calculated in the previous layers. Analysis towards the center leads the evolution of the reconstruction. This multi-layer approach is adapted to the tree because their circumferential cross section, from internal to external, can be associated to a 2D annular structure.

Discretization

The number of layers depends on the number of transducers. N (pair) transducers are distributed uniformly around the object every 360 o / N (angular step). Equation 5shows that the value of the radius of a layer is determined by:

d(k) = min M∈D |C -M| ; ∀ 1 ≤ k ≤ N/2; k ∈ N (5)
Where: d(k): the distance from the layer k towards the center (radius of the layer). Thus the mesh is composed by N/2 layers. C: is the center of the object. D: the straight ray linking transducers E 1 E k+1 M: point in D

Algorithm

Let Eq. 6 describe the slowness for a simple backprojection process. Γ is the straight ray connecting the emitter E(n) to the receiver E(m). Ω(i, j) the set of straight rays that intersect a pixel at (i,j), where i and j are the index of the resulting image by the reconstruction 1 . |Ω| is the number of straight rays of the set.

p(i, j) = 1 |Ω| Ω(i, j) P(n, m) (6)
From the foward problem, and assuming the existance of nested layers, Eq. 4 could be explained by Eq. 7 .

NP(n, m) = layer k-1 p k-1 (i, j) + layer k p k (i, j) (7)
A generalization of Eq. 4 leads to the slowness P of a layer. Slowness at a layer k-1 and k is then estimated as shown in Eq. 8 and Eq. 9.

layer k-1 p k-1 (i, j) = N k-1 P k-1 ( 8
)
layer k p k (i, j) = N k P k (9)
where N k-1 and N k are the number of pixels at layers k-1 and k respectively.

Replacing Eq. 9 in Eq. 7, the correction formula for P(n,m) is then calculated as follows:

NP(n, m) = layer k-1 p k-1 (i, j) + N k P k (10) Pk (n, m) = NP(n, m) - layer k-1 p k-1 (i, j) N k (11)
where Pk (n, m) is the corrected value of the backprojection. Recursive application of Eq. 11 into Eq. 6 leads to the reconstruction. For each iteration, the upper part of Eq. 11 should be positive to guarantee the convergence to a solution.

Results

Forward Problem

Figure 3 shows the ultrasonic sinogram formed by 90 observation angles (angular step of 4 o ). Ray theory was applied without taking into account the attenuation in the different media. The color of each pixel corresponds to the slowness measured between two elementary sensors. Each line of the matrix corresponds to a projection.

Reconstruction and processing -Inverse Problem

Figure 4 shows the reconstructed object associated to each inversion algorithm. The reconstructed objects are smaller with respect to the original phantom. The theoretical number 1 Size of the phantom and reconstructed image are different.

of pixels that allow to reconstruct the scanned area with 90 emitter positions is 4005 [START_REF] Socco | Feasibility of ultrasonic tomography for non-destructive testing of decay on living trees[END_REF]. The optimal size of the associated image is then 72x72 pixels. Note that the quality of the reconstructions obtained in Figure 4 is partly due to the fact that the forward problem and inverse problem involve the same approximate theory of straight ray propagation. Figure 5 shows the uncertainty of the linear model Y=aX+b that compares the reconstructed measured slowness (X = projections of the reconstructed image (depend of the algorithm)) with respect to initial measured values (Y = projections of the phantom). The analysis were carried out with a confidence interval of 95% and each test was realized for the reconstruction for a given number of sensors. Figure 6 shows the mean of absolute deviations (mean deviation per pixel) between the reconstructed image and the phantom. The analysis was carried out by linear interpolation of the reconstructed image to the size of the phantom and comparing the deviation pixel per pixel. When noise is not taken into account, both figures show that when this range of sensors (24 -96) is used, the uncertainties decrease. Because the observations for further than 96 sensors were not tested, conclusion about the asymptotic behavior was not established. Reconstructions by the Layer Stripping algorithm are less performing than those carried out by the inversion by FBP. For a scan with 90 sensors, the respective uncertainties of the linear model were 7.2% (FBP) and 12.6% (LS). The mean deviations per pixel were about 0.079x10 -3 sm -1 (FBP) and 0.124x10 -3 sm -1 (LS); to be compared with the mean pixel value of 0.92x10 -3 sm -1 inside the scanned area of the phantom. For each one of these figures the slopes of each pair of curves are similar. Both inversion methods present similar behavior according to the number of sensors used for a given scan. 7 and8 show respectively the evolution of the linear model uncertainty and of the mean deviation per pixel when Gaussian noise is present for the same scan with 90 sensors. As noise increases, the uncertainties do it as well, but the behavior for the two inversion methods is different. The slopes of the curves of the uncertainties for the reconstructions by FBP are more sensitive to noise than for those by Layer Stripping. From about 20% of noise, the mean deviation per pixel for Layer Stripping is lower than for the FBP. For example at 30% of noise, linear model uncertainties are 15.7% (FBP) and 16.8% (LS), and are related to 0.209x10 -3 sm -1 and 0.173x10 -3 sm -1 deviation per pixel respectively.

Conclusion

A new inversion algorithm has been presented as well as first numerical results. The tests consisted to simulate travel time tomography with the aim to compare the inversion of numerical data by Layer Stripping in relation to the inversion by Filtered Backprojection (reference method). These simulations depended of the experimental set-up of the ultrasonic tomograph and assumptions for the numerical phantom. For this first approach the phantom was assumed to be a cylindrical homogeneous isotropic medium with defects. Ray theory was applied without taking into account the attenuation. Quantitative analysis included the comparison of the uncertainties of the projections and slowness values of the reconstructed images with respect to the phantom. The first analysis showed that if noise is not taken into account, as more sensors are used for a scan, uncertainties decrease. Both inversion algorithms presented similar behavior according to the number of sensors. A general conclusion about this behavior was not established because tests were developed for a specific range of number of sensors. The results showed that if noise is not taken into account, uncertainties for Layer Stripping are greater than for FBP. The second analysis showed that, for a given scan, as noise increases, uncertainties do it as well, however the reconstructions by the FBP became more sensitive than those by Layer Stripping.

These results are the first step for further research in numerical modeling and experimental testing with samples of trees. In this paper it was underlined that the approach by Layer Stripping could be adapted to trees because it can be associated to their 2D annular structure. A next step in numerical modeling should be the choice of a phantom with a concentric layer geometry. In relation to the inversion algorithms, because the actual challenge in tree imaging is the choice of the effective method to inverse scattered data, several reconstruction methods, from analytical to algebraic, should be also tested. Experimental tests using our tomographic device are also planned to be carried out.

Figure 2 :

 2 Figure 2: Numerical phantom. o : Positions where the sensors are going to be experimentally placed.

Figure 3 :

 3 Figure 3: Simulated ultrasonic sinogram.

Figure 4 :

 4 Figure 4: Reconstruction by: a) Filtered backprojection; b) Layer Stripping

Figure 5 :

 5 Figure 5: Linear model uncertainty according to the number of sensors for: o Filtered backprojection, * Layer Stripping

Figure 6 :

 6 Figure 6: Mean deviation per pixel according to the number of sensors for: o Filtered backprojection, * Layer Stripping

Figures

  Figures7 and 8show respectively the evolution of the linear model uncertainty and of the mean deviation per pixel when Gaussian noise is present for the same scan with 90 sensors. As noise increases, the uncertainties do it as well, but the behavior for the two inversion methods is different. The slopes of the curves of the uncertainties for the reconstructions by FBP are more sensitive to noise than for those by Layer Stripping. From about 20% of noise, the mean deviation per pixel for Layer Stripping is lower than for the FBP. For example at 30% of noise, linear model uncertainties are 15.7% (FBP) and 16.8% (LS), and are related to 0.209x10 -3 sm -1 and 0.173x10 -3 sm -1 deviation per pixel respectively.

Figure 7 :Figure 8 :

 78 Figure 7: Linear model uncertainty according to the percentage of noise considered for: o Filtered backprojection, * Layer Stripping

Table 1 :

 1 Parameters of the phantom

	Diameter (in pixels)	Position	Slowness Velocity (10 -3 sm -1 ) (ms -1 )
	122	centered	0.74	1350
	28	off-centered	1.11	900
	58	off-centered	1.43	700

Acknowledgments

This study was based on the research supported by the French National Research Agency (BioGMID -Biological Growth Medium Integrity Diagnoses using bi-modality tomographies) under the Grant n o 183692 to USAR-CNRS.


For cortical bone, important changes of the elastic properties values have been clearly shown in ageing but not in childhood, furthermore recent works considered osteoporosis as a pediatric disease with geriatric consequences and children are concerned by specific infantile osteo-pathologies. That is why there is a strong interest in the characterisation of the growing process of children bone. However, few mechanical properties of cortical growing bone are available in literature and do not yield to gold standards. Results found in litterature for children bone concern specific location (close to cancerous cells) or cadaveric bone. They indicate a lower Young's modulus for children bone compare to mature bone. The goal of this study is to provide elastic properties values for human growing cortical bone. To reach this goal, we have analysed surgery waste (bone transplantation) from long bone (fibula). In a first step, a non destructive method was used to evaluate the velocity of ultrasonic waves from which the acoustic Young's modulus E a is calculated using the difference of sound path duration and the mass density. Then, in a second step, a destructive method was used to obtain mechanical Young's modulus E m using a 3-point microbending. Heihteen children bone samples were tested acoustically and due to the specific dimensions needed for micro three points bending fourteen samples were tested this way, but for each step, comparisons were made with adult specimens (twelve samples). The children surgery wastes from fibula (4 to 16 year old children) included in this study show an average E a and ν a of 15.5 GPa (+/-3.4) and 0.24 (+/-0.08) at 10 MHz , and an average E m of 9.1 GPa (+/-3.5). E a and ν a are in the same range for children and seniors but a linear correlation between E a and E m is found only for the fourteen samples of the children group.

Introduction

The cortical bone constitutes 80% of the human skeleton mass, hence its characterisation represents a main stake particularly for the diagnosis of several pathologies and most of the time on aging pathologies. Indeed, aging adversely affects the elastic properties of human cortical bone throughout life, from childhood to old age. Lot of studies highlighted differences in the mechanical behavior of physiological and osteoporotic adult bone but the evolution in mechanical properties of bone in childhood has been less studied and nowadays there is a tremendous lack of data on human growing bone.

In paediatrics, diagnosis of bone diseases and treatment decision are provided without specific knowledge of bone mechanical properties. The paediatric numerical models used in orthopaedics or in accidentology include adult mechanical parameters or values from a theorical method displayed by [START_REF] Irwin | Biomechanical basis for the crabi and hybrid iii child dummies[END_REF]. This latter is based on the hypothesis of a decrease in Young's modulus from adult values to neonate values. But, this theorical approach have been invalidated for the Young's modulus of the skull bone [START_REF] Roth | Limitation of scaling methods in child head finite element modelling[END_REF].

To assess the children bone quality, the quantitative ultrasound techniques seem promising as they are non-radiating and non-invasive. Nevertheless, in acoustical imaging, there is no data for children cortical bone to provide quantitative information about the mechanical characteristics of the children tissues analysed. For cortical bone, few studies considered mechanical characteristics of the growth process [START_REF] Currey | The effect of porosity and mineral content on the Young's modulus of elasticity of compact bone[END_REF]. Moreover in most of the case they have chosen specific bone, close to cancellers cells, or cadaver fragments [START_REF] Baleani | Compressive Behaviour of Cortical Bone in Young Humans[END_REF][START_REF] Chotel | Bone stiffness in children: part i. in vivo assessment of the stiffness of femur and tibia in children[END_REF], and they did not provide mechanical information on physiological bone.

The aim of this study is to give values of elastic properties (E and ν) of the human growing bone.

In this study, we have used surgery wastes (bone transplantation) coming from fibula (long bone) essentially composed of cortical bone. A two step method has been carried out, first an ultrasonic evaluation which provided the transverse Young's modulus and Poisson's ratio, and second, a three point microbending test which provided longitudinal Young's modulus. These values have been compared to a senior group composed of mature bone fibula (+ 75 YO).

Materials and Methods

Samples

Using a Marseilles hospital population, we studied cortical bone from autotransplants. From non pathological location, the whole autotransplant is extracted at the bottom of the fibula 5 cm up to the ankle. The non-use parts (wastes) are kept and cut using a low speed diamond saw (Isomet 1000, Buehler; Lake Bluff, IL, USA) in order to obtain parallelepipedic samples (plane and parallel surfaces). A total of 30 rectangular samples with parallel edges were cut. Their dimensions were measured with a digital calliper. One of the major difficulty looking at children bone is the small size of the samples: 15 to 35 mm long (bone axis direction), 10 to 20 mm wide and 2 to 3.5 mm thick (transverse directions). They are designated as (F or M for sex)-(age)-(F for fibula)-(number for each piece) and stored at -20 • C. Two groups have been studied : 18 samples composed the children group (4-16 year old) and 12 samples constitute the senior group (+75 year old).

First step : ultrasound measurements

A specific ultrasound frame has been developed to process very small samples. The ultrasonic experimental setup is presented on Figure 1. It provides the evaluation of lon-10 Figure 1: Ultrasonic experimental setup gitudinal and transverse waves velocities in transverse direction (perpendicular to the bone axis) from the measurement of the time-of-flight (TOF) before and after the first critical angle. This protocol has been detailed in [START_REF] Bensamoun | Transmission scanning acoustic imaging of human cortical bone and relation with the microstructure[END_REF] The samples prepared for ultrasonic measurements were completely immersed in a temperature-controlled tank filled with distilled, degassed water at 25 • C. They are mounted with y-axis corresponding to the axis of the bone (longitudinal direction), and x-axis and z-axis are the radial and circonferential directions (transverse directions). They are moved with 0.1 mm accuracy in x and y-axis. For z-axis, azimuth got 0.01 • of accuracy. Two focused broadband transducers at 10 MHz (5 mm diameter, focal length 30 mm ; Imasonic, Besanc ¸on, France) operated through each sample. The sample thickness was first calculated using the echo-technique with transducers used as both the transmitter and receiver and then compared to calliper measurements.

Then, the TOF of longitudinal and transverse waves was measured in transmitter-receiver mode on a mechanical support allowing translations in 3D. Given the thickness of the sample and the distance between the transducers, TOF measurement at normal incidence gives the longitudinal wave velocity (V L ) and the TOF measurement after the first critical angle gives transverse wave velocity (V T ) as done by [START_REF] Pithioux | An alternative ultrasonic method for measuring the elastic properties of cortical bone[END_REF].

Finally, using the Archimedes' principle, a micrometric balance (610 GX0.001G voyager) measured bone mass density, ρ (each parameter was assessed by performing tested triplicate). Therefore, the transverse Young's modulus E a and Poisson's ratio ν a were provided through the following equations:

E a = ρ * V 2 T * (3V 2 L -4V 2 T ) (V 2 L -V 2 T ) (1)
ν a = V 2 L -2V 2 T 2(V 2 L -V 2 T ) (2)

Second step : 3-point microbending testing

A specific 3-point microbending testing system has been designed for small samples and mounted on an Universal Testing Machine(Instron 5566A, Norwood, MA). To evaluate the cortical bone samples, a span to depth ratio of 10:1 was respected which reduced the number of tested samples to 14 for the children group. The smallest length was 15 mm and the smallest thickness 1.5 mm. A pre-force of 5 N was applied on the sample before testing till rupture. The displacement speed was 0.1 mm/min, close to static testing conditions; the test provided a force/displacement curve for each sample which is transformed in strain/stress curve from which the Young's modulus, noted E m , is estimated. age groups (Table 1). The measurement error of ultrasound velocity is estimated to 2.25%.

Results

Ultrasound measurements

Three-point microbending test

Three-point microbending test provides an average longitudinal Young's modulus E m of 9.1 GPa (+/-3.5) for the children group and of 5.8 GPa (+/-2.1) for the senior group (Table 1

Discussion

Both the values of E a and E m are consistent with values usually reported in the litterature [START_REF] Dong | The dependence of transversely isotropic elasticity of human femoral cortical bone on porosity[END_REF][START_REF] Grimal | Assessment of cortical bone elasticity and strength: Mechanical testing and ultrasound provide complementary data[END_REF].

The elastic properties of children bone assessed by ultrasonic method are compared to old age bone samples extracted from the same anatomical site (bottom of the fibula, 5 cm up to the ankle). Our results show no statistic difference in transverse Young's modulus and Poisson's ratio between children and old age bone (Student T test p>>0.05). The adult bone values found in this study spread in the typical range found in the literature [START_REF] Yoon | Ultrasonic wave propagation in human cortical bone -ii measurements of elastic properties and microhardness[END_REF][START_REF] Berteau | Caractérisation mécanique de l'os en croissance par méthode ultrasonore[END_REF]. Concerning the Poisson's ratio, the average values is roughly 0.2-0.3, as usual for cortical bone. Furthermore, mass density of cortical bone from fibula measured according to the Archimedes' principle (1)are in the usual range found in the literature [START_REF] Zioupos | Some basic relationships between density values in cancellous and cortical bone[END_REF]. The values obtained for the mechanical properties are quite similar between children bone samples and old age bone samples. Nevertheless, further experiments are needed to explore the age-related variation in mechanical properties, in particular the age ranges between the children group and the senior group. Indeed, refering to the paper of [START_REF] Drozdzowska | Skeletal status in males aged 7-80 years assessed by quantitative ultrasound at the hand phalanges[END_REF], the speed of sound (SOS) evaluated at the phalanx is quite the same for children around 10 year old and 70-80 year old people. But Drozdzowska and colleagues put on evidence that the SOS at the phalanx increases linearly to a maximum value reached around 25 year old and then the values decrease more slowly up to the 80 year old.

These acoustical measurements on mechanical properties of children bone are looked at in relation to the results obtained from three-point microbending test (Figure 3). It has to be noted that the span to depth ratio of 10:1 is not theoretically sufficient and can induce some errors in the evaluation of E m as shown in [START_REF] Spatz | Young's Moduli and Shear Moduli in Cortical Bone[END_REF]. Nevertheless, all the samples are in the same span to depth ratio which allows to compare the results obtained for the children group to the results obtained for the senior group. We investigate the correlation between the elasticity in the two directions: longitudinal (E m ) and transverse (E a ). In the literature, some studies on mature bone found low positive or non-signifant correlations between elasticity in longitudinal and transverse directions. A recent study (Gri-mal2009) obtained a marginal negative correlation in agreement with the results we obtained on the senior group (Figure 3 (b)). On the other hand a very different trend appears concerning the 14 children bone samples testing by the two modalities: Figure 3 (a) shows a high positive corelation between E a and E m . To our knowledge this result is the first one showing a strong link between the elasticity of the longitudinal and transverse directions for human growing bone. Further work is needed to extend these observations to the anisotropy of growing bone which is a major parameter of the mechanical behaviour of children bone but still unknown.

Conclusion

Even if changes of material properties of cortical bone are more and more studied in ageing, this changes are not clearly shown in childhood. Indeed, some studies compare bone density (Bone Mineral Density : BMD) of children, like in Crohn's disease [START_REF] Semeao | Bone mineral density in children and young adults with Crohn's disease[END_REF], but not the stiffness, strength or elastic modulus. The few studies considering mechanical characteristics of this process have chosen specific bone, close to cancerous cells [START_REF] Ohman | Compressive behaviour of child and adult cortical bone[END_REF], or cadaver fragments [START_REF] Currey | The mechanical properties of bone tissue in children[END_REF]. One of the goal of this study is to underline that the use of mature bone elastic properties to study growing bone pathologies is unsuitable. As a consequence the work started in this study is particularly relevant to gain insights on the specific mechanical behaviour of children bone.

Figure 2 Figure 2 :

 22 Figure 2 reports the values of longitudinal and transverse wave velocities calculated for the 18 children bone samples. The average values of longitudinal velocities are 3225 m/s for the children group and 3499 m/s for the senior group and the average values of transverse velocities are 1827 m/s for the children group and 1930 m/s for the senior group.The deduced transverse Young's modulus and Poisson's ratio, E a and ν a , are around 17 GPa and 0.25 respectively, for the two

Figure 3 :

 3 Figure 3: Correlation between Young's modulus measured by 3-point bending test (longitudinal direction) and Young's modulus measured by ultrasound (transverse direction): a)children group; b)senior group

Table 1 :

 1 ). The measurement error of the cell-force is estimated to 0.23.% Average values and standard deviation of E a , ν a and E m comparing children group samples to senior group samples

		ρ	E a	ν a	E m
		kg/m 3	(GPa)		(GPa
	Children	1.84	15.5 ± 3.4 0.24 ± 0.08 9.1 ± 3.5
	Senior	1.73	16.7 ± 1.9 0.27 ± 0.05 5.8 ± 2.1

Acknowledgments

We are grateful for financial assistance to the ANR project called Biogmid, we thank the Timone's Hospital surgery team and the donors or their legal guardians who gave informed written consent to providing their tissues for investigation, in accordance with the French Code of Public Health (Code de la Santé Publique Franc ¸ais) and approved by the Committee for the Protection of Persons.


Hamlet Medina Ruiz

Michel Kieffer

Béatrice Pesquet-Popescu

An Adaptive Redundancy Scheme for TCP with Network Coding

To address the inability of the standard TCP protocol to distinguish between congestive losses and random packet losses produced by the noisy channel, in this paper we propose an adaptive algorithm to dynamically adjust the redundancy factor R of the TPC/NC protocol proposed by Sundararajan et al. by adding some additional functionalities to the original network coding layer. We define a loss differentiation scheme to adjust R, based on the Vegas Loss Predictor and the collective feedback information of ACKs and duplicates ACKs, which are both indicators of the network condition. In that way the source adjusts intelligently R based on the network conditions, avoiding unnecessary TCP rate reduction due to packet losses induced by the noisy channel and preventing the network from entering in a congestion state. Our TPC/NC implementation with our adaptive scheme is full-duplex, and can manage multiple TCP connections simultaneously. Simulation results over realistic network scenarios using

OPNET Modeler show that our adaptive scheme in conjunction with the standard TCP/NC produces better TCP-throughputs than the standard TCP/NC, TCP-Reno, TPC New Reno, and TCP Reno with SACKS.

I. INTRODUCTION

Network Coding (NC) is a technique that offers a lot of potential in today's applications by allowing nodes to send out packets that are linear combinations of previously received information, instead of delivering the information to their destination in the standard store-andforward-manner [START_REF] Ahlswede | Network information flow[END_REF], [START_REF] Robert Li | Linear network coding[END_REF]. The main benefits of NC are the potential throughput improvements and a high degree of robustness, which is translated into loss resilience [START_REF] Fragouli | Network coding: an instant primer[END_REF]. Despite all the benefits of NC , we still lack of many real network coding implementations.

In [START_REF] Sundararajan | Network coding meets TCP[END_REF] Sundararajan et al. propose a new TCP-friendly protocol that interfaces network coding with TCP by introducing a new coding layer in the TCP/IP stack between TCP and IP, with the idea of using the benefits of NC and to improve TCP throughput. In this layer TCP segments are encoded at the sender and decoded at the receiver. In [START_REF] Sundararajan | Interfacing network coding with TCP: an implementation[END_REF] the authors also describe several important practical aspects of incorporating NC with TCP-Reno, as the TCP-compatible sliding window code, and new rules for acknowledging bytes to the TCP layer. Redundancy of information added by the sender is used to cover non congestion losses from TCP congestion control, and simulation results in [START_REF] Sundararajan | Network coding meets TCP[END_REF], [START_REF] Sundararajan | Interfacing network coding with TCP: an implementation[END_REF] show an improvement in TCP throughput for a constant value for R.

Setting R without having an estimation of the network conditions, as the packet loss rate of the channel and the congestion in the network at any point in time, could lead to TCP times outs when R is too low to match TPC's sending rate to the rate at which data is received at the receiver because of the losses, or, on the other hand, if R is too large, the losses can be [START_REF] Floyd | The NewReno modification to TCP's fast recovery algorithm[END_REF] DRAFT recovered, the congestion in the network could increase and the TCP throughput decrease . As the correct approach is to increase the transmission rate when the channel is dropping packets [START_REF] Tanenbaum | Computer Networks[END_REF], [START_REF] Peterson | Computer Networks: A Systems Approach[END_REF], we would like to increase R in these situations to compensate losses, as for example when TCP-Reno probing for available bandwidth believes that it has reached the network congestion state when it wrongly identifies the reception of duplicates ACKs (DUPACKS) caused by the noisy channel as a sign of congestion, resulting in an unnecessary reduction of the transmission rate, and in a similar manner decrease R when the reception of DUPACKS are really caused by a congested network.

We propose a new adaptive algorithm to adjust R that uses the feedback information available at hand, namely packet drops and packets backlogged based on the Vegas Loss Differentiation algorithm [START_REF] Brakmo | TCP vegas: end to end congestion avoidance on a global Internet[END_REF] to define a loss differentiation scheme, and based on this information adjust R, helping the TCP congestion control mechanism to properly reduce its congestion windows only when there exists congestion in the network, and reacts to packets drops in a smooth manner.

This paper evaluates the performance of TCP/NC with adaptive R vs the standard TCP/NC protocol, as well as with the standard TCP congestion control flavors, namely TCP-Reno, arrived packets in the coding buffer. To convey information about the linear combination, the NC header is added to the coded packet, which contains information about the coefficients used to mix the packets, the starting and last byte SEQ of each packet involved in the mix, and the first byte that has not been ACKed. The original packets are retained in the coding buffer until an appropriate TCP ACK arrives from the receiver side.

On the receiver side, upon receiving a new linear combination, the decoder places it in a decoding buffer, appends the corresponding coefficients to the decoding matrix ( which contains the coefficients used to mix already seen packets received), and performs Gaussian elimination.

This process helps identify the "newly seen" packet (if it exists). Then the receiver sends a TCP cumulative ACK to the sender requesting the first unseen packet in order, like in conventional TCP. The Gaussian elimination may result in a new packet being decoded. In this case the decoder delivers this packet to the receiver TCP. Any ACK generated by the receiver TCP is suppressed and not sent to the sender. These ACKs may be used for managing the decoding buffer. An important point is that the introduction of the new NC layer does not require any change in the basic features of TCP. For more detail, the reader is referred to [START_REF] Sundararajan | Interfacing network coding with TCP: an implementation[END_REF].

III. ADAPTIVE TCP/NC PROTOCOL

A. Adaptive redundancy factor

The main idea behind coding is to mask the losses on the channel from TCP. In other words, we wish to correct the losses without relying on the ACK's.

Linear combinations generated using NC should be able to correct any loss of the W packets involved in the linear combinations that have not yet been ACKed. R linear combinations on average are sent to the IP layer to compensate for the loss rate of the channel, and to match TCP's sending rate to the rate at which data are actually sent to the receiver. If there is too little March 9, 2012 DRAFT redundancy, we will not receive enough linear combinations to decode (losses are not masked from TCP), leading to timeouts and consequently low throughput.

On the other hand, if the redundancy is too large, losses are recovered, but TCP congestion window advances smoothly, the throughput is reduced due to low code rate and the congestion increases.

Assuming a constant redundancy factor equal to the reciprocal of the probability of successful reception is impractical, as this probability changes over time due to changes on bit error probability of the channel and the congestion in the network.

The main idea of this paper is to adjust dynamically R to an optimal rate based on the actual network conditions, which can be detected by collecting the feedback information available. Our objective is to adjust R in such a way that TCP congestion window slows down its sending rate vigorously to alleviate congestion only when we have an indication of a really congested network, and in the worst case to maintain its rate when loss packets are due to the noisy channel . When the channel is dropping packets we try harder by increasing the redundancy to mask the channel losses from TCP, and we decrease the redundancy when congestion is present, allowing TCP to sense the congestion and reduce its rate to alleviate it. To accomplish this objective, we make use of the TCP Vegas loss predictor to detect when congestion is present in the network, in conjunction with the number of duplicate ACK's generated by the receiver.

1) Vegas loss predictor: TCP Vegas uses a loss predictor to decide whether the network is congested or not based on rate estimators [START_REF] Brakmo | TCP vegas: end to end congestion avoidance on a global Internet[END_REF], [START_REF] Martignon | Loss differentiation schemes for TCP over wireless networks[END_REF], [START_REF] Tian | TCP in wireless environments: problems and solutions[END_REF] , estimating the backlogged packets in the buffer of the bottleneck link. We implement the Vegas Loss Predictor at the NC layer to know when the network experiences congestion and to adjust R in a more intelligent way. As in [START_REF] Brakmo | TCP vegas: end to end congestion avoidance on a global Internet[END_REF], we set BaseRT T to be the RTT of a segment when the connection is not congested (in March 9, 2012 DRAFT practice the minimum of all measured round trip times, which is reset after a time out) as a reference to derive the expected throughput the network can accommodate.

The expected throughput is given by window/BaseRT T , where window is the difference between the SEQs of the newest and oldest byte involved in the linear combinations currently in transit.

The measured throughput is given by windowAcked/RT T , where windowAcked is the number of bytes ACKed (seen by the receiver NC layer) during the last RTT, and RT T is the sample RT T of the received segment.

To estimate the cause of the packet losses, the parameter qv is calculated as

qv = window BaseRT T - windowAcked RT T BaseRT T (1)
and compared with two thresholds α and β to know the network state. If qv ≥ β, the network is considered congested (too much backlogged), if qv ≤ α, possible losses due to random errors (little backlogged), and finally if α < qv < β, the predictor assumes that the network state is the same as in the previous estimation. In our test the best performance was achieved with the common values used in Unix OS flavor implementations, α = 1 and β = 3.

A detailed accuracy analysis of the Vegas predictor and other loss differentiation algorithms can be found in [START_REF] Martignon | Loss differentiation schemes for TCP over wireless networks[END_REF].

2) New Adaptive algorithm for R: In this algorithm we define a loss differentiation scheme with the help of the Vegas Loss Predictor [START_REF] Brakmo | TCP vegas: end to end congestion avoidance on a global Internet[END_REF], and the number of DUPACKS generated by the receiver to detect if a packet loss is due to the noisy channel or to a congested network by comparing these two measures.

We propose an implicit approach to adjust R based on this differentiation scheme, allowing the sender to react intelligently to the cause of the packet loss, due to either congestion or random March 9, 2012 DRAFT errors, and therefore prevents the network from entering in congestion or an unnecessary decrease of the congestion window.

Basically, each time an ACK is received, the quantity in (1) is evaluated. If the value of ( 1) is below the threshold that indicates a congested network and some DUPACKS have been received, it is because packets are being dropped by the noisy channel and not by buffer overflows in the router queues'. In these situations the correct approach is to increase R to mask the losses. On the other hand, if a network congestion is detected using (1), we decrease R, allowing TCP to also reduce its rate and alleviate the congestion. Initially we set R to a value that takes into account the loss in throughput due to the finiteness of the field to reduce the number of DUPACKS due to dependent linear combinations received at the receiver, and when congestion is detected R is set again to this value. The algorithm is specified below using pseudo-code.

In order to successfully implement Algorithm 1 in the NC layer, some minor improvements are added to the standard TCP/NC protocol. Timer functions regarding the RTT calculation, and variables to keep track of the number of loss events based on the reception of DUPACKS, the value of window and windowAcked in (1) have been implemented.

Our implementation of the TCP/NC protocol with the proposed algorithm is full-duplex, and also supports multiple instances, one for each TCP connection present [START_REF] Stevens | TCP/IP illustrated (vol. 2): the implementation[END_REF], which is identified by the TCP socket. To put it clear, the algorithms for the source and receiver side are shown is pseudo code below. This specification assumes only one instance and one-way TCP flow. The source side algorithm has essentially to respond to two types of events, the arrival of a packet from the source TCP, and the arrival of an ACK from the receiver via IP. On the receiver side, the algorithm has to respond to the arrival of a packet from the source, and the arrival of ACKs from the TCP sink.

if qv ≥ β then R ← 1.1 else if qv ≥ α then if cntDupAck ≥ 2 then R ← R + 2/10; else if cntDupAck == 1 then R ← R + 1/10 else R ← R -1/10 end if else R ← 1.1 end if

IV. SIMULATION SCENARIO AND RESULTS

A. Simulated Topology

We use a client/server OPNET model to evaluate and compare the performance of TCP/NC with adaptive redundancy vs the standard TCP/NC protocol, TCP-Reno, TCP-NewReno and TCP-Reno with SACK congestion control mechanisms [START_REF] Jacobson | Congestion avoidance and control[END_REF], [START_REF] Fall | Simulation-based comparisons of tahoe, reno and sack TCP[END_REF], [START_REF] Floyd | The NewReno modification to TCP's fast recovery algorithm[END_REF], in a scenario when channel and congestion losses are simultaneously present.

B. Simulation Setup

In this scenario the network carries traffic from three different applications, FTP, Video, and VoIP. VoIP traffic has PCM quality speech, and Video traffic is a low resolution video streaming, with a frame inter-arrival time of 10 fps and frame size of 128 × 120 pixels.

Each type of traffic is marked with a priority using the IP Type of Service (ToS) field. In this scenario the FTP application begins the connection at time t = 10 s and finishes when the file is totally transferred or by the end of the simulation time which is of 15 minutes. The VoIP and Video Streaming traffic are active between 60 s < t < 110 s and between 80 s < t < 110 s, respectively, generating congestion in the network during that period of time. In the simulation we have also modeled the loss in throughput due to the finiteness of the field (in our case, its size is 256).

March 9, 2012 DRAFT We assume for the rest of the paper a simulation time of 15 minutes, a latency equal to 100ms on the DS1 link, a FIFO queueing discipline with a buffer space equal to 500 in each router interface queue (value that is much larger than the bandwidth-delay product) , and the TCP simulation parameters shown in Table I. The parameters in Table I are used for all the TCP flavors.

Performance parameters as the congestion window size, file download response time and goodput are studied through the simulations.

March 9, 2012 DRAFT In Table II we also compare the download time response experienced by the FTP client since the FTP connection is established, in the case when the loss rate is equal to 3%.

In Figure 6 Finally, in Figure 7 we show the evolution of the TCP congestion window (cwnd in bytes) as a function of the redundancy factor R during first 450s of simulations for a loss rate of 3%, using adaptive R and contant R equal to 1.3. We observe that at the moments when cwnd curve presents almost flat segments there are peaks in redundancy curve. This represents events when DUPACKS are received, caused either by the noisy channel or by the finiteness of the field. It is also important to notice, that our scheme adapts well to a burst of errors, which can be viewed by comparing the oscillating behavior of the cwnd for R = 1.3 and the smoothness of the cwnd using our scheme in some time intervals, without having a too large value for R. Increasing R when DUPACKS are received by non congestion causes prevents the unnecessary reduction of TCP sending rate.

V. CONCLUSIONS AND PERSPECTIVES

This paper considers the issue of adapting over time the redundancy factor R in the original TCP/NC protocol.

We proposed a new adaptive algorithm to dynamically adjust R using packet drops and packets backlogged information based on received DUPACKS and the Vegas Loss Differentiation algorithm, respectively, to discover the underline cause of a packet loss in the network and take the correct action over R.

Simulation results have shown an improvement in the goodput and download FTP time In future work, simulations over heterogenous networks, a mathematical proof of the stability of the TCP/NC protocol with our adaptive redundancy scheme, as well as a characterization of the throughput behavior as a function of the loss rate will be considered.

March 9, 2012 DRAFT

Algorithm 1

 1 Calculate qv α = 1, β = 3, R = 1.1 {Initialization}Each time an ACK is received: Call timer routines and compute qv[START_REF] Ahlswede | Network information flow[END_REF]

3Fig. 1 .

 1 Fig. 1. Topology of the wide Area Network

FTP traffic is configuredFig. 2 .

 2 Fig. 2. Servers and VoIP client located in subnetA

Fig. 3 .

 3 Fig. 3. Clients Located in Subnet B

Fig. 4 .

 4 Fig. 4. Goodput vs Loss Rate for TCP/NC and TCP/NC with adaptive R

Fig. 5 .

 5 Fig. 5. Goodput vs Loss Rate for TCP/NC with adaptive R and different TCP Flavors

Fig. 6 .

 6 Fig. 6. Sent SEQ vs Time for TCP/NC and TCP/NC with adaptive R

Fig. 7 .

 7 Fig. 7. Evolution of the redundancy and the congestion window for 1 % loss probability

ACKNOWLEDGMENTS

This work has been partly supported by DIM-LSC NC2, DIM-LSC SWAN. Michel Kieffer was partly supported by the Institut Universitaire de France.

Algorithm 3 Receiver Side Algorithm

1 Wait state: If any of the following events occurs, respond as follows; else, wait.

2 ACK arrives from TCP sink: If the ACK is a control packet for connection management, deliver it to the IP layer and return to the wait state; else, ignore the ACK.

C. Simulation Results

Using the simulation setup explained in section IV-B, we first study the effect of the parameter R on the goodput of the standard TCP/NC protocol and our TCP/NC approach as a function of the loss rate of links (the loss rate being the probability of a packet getting lost on the link). Packets in the network are subject to these losses in the forward and the reverse direction.

Figure 4 shows the plot of the TCP/NC goodput as a function of the loss rate when the parameter R is varied between 1.0 and 1.5, as well as when R is dynamic. From this plot, it is clear the important role played by parameter R. We can see that for the range of R and loss rates tested, when R is varied over time according to Algorithm 1, we outperform the standard TCP/NC .

In order to compare the performance obtained with our approach with the standard TCP congestion mechanisms, in Figure 5 is plotted the TCP goodput for TCP-Reno, TCP-NewReno, March 9, 2012 DRAFT