A parametric imaging approach for the segmentation of ultrasound data

Ultrasonics 43 (2005) 789–801 www.elsevier.com/locate/ultras

A parametric imaging approach for the segmentation of ultrasound data F. Davignon *, J.-F. Deprez, O. Basset Creatis UMR CNRS 5515, Unite´ INSERM 630, Univ. C. Bernard, Baˆtiment Blaise Pascal, 69621 Villeurbanne Cedex, France Received 28 January 2005; received in revised form 5 June 2005; accepted 12 June 2005 Available online 6 July 2005

Abstract When an ultrasonic examination is performed, a segmentation tool would often be very useful, either for the detection of pathologies, the early diagnosis of cancer or the follow-up of the lesions. Such a tool must be both reliable and accurate. However, because of the relatively reduced quality of ultrasound images due to the speckled texture, the segmentation of ultrasound data is a difficult task. We have previously proposed to tackle the problem using a multiresolution Bayesian region-based algorithm. For computation time purposes, a multiresolution version has been implemented. In order to improve the quality of the segmentation, we propose to perform the segmentation not only from the envelope image but to combine more information about the properties of the tissues in the segmentation process. Several acoustical parameters have thus been computed, either directly from the images or from the radiofrequency (RF) signal. In a previous study, two parametric images were involved in the segmentation process. The parameter represented the integrated backscatter (IBS) and the mean central frequency (MCF), which is a measurement related to the attenuation of ultrasound waves in the media. In this study, parameters representative of the scattering conditions in the tissue are evaluated in the multiparametric segmentation process. They are extracted from the K-distribution (a,b) and the Nakagami distribution (m,X) and are related to the local density of scatterers (a,m), the size of the scatterers (b) and the backscattering properties of the medium (X). The acoustical features are calculated locally on a sliding window. This procedure allows to built parametric mapping representing the particular characteristics of the medium. To test the influence of the acoustical parameters in the segmentation process, a set of numerical phantoms has been computed using the Field software developed by J.A. Jensen. Each phantom consists in two regions with two different acoustical properties: the density of scatterers and the scattering amplitude. From both the simulated RF signals and envelope images, the parameters have been computed; their relevance to represent a particular characteristic of the medium is evaluated. The segmentation has been processed for each phantom. The ability of each parameter to improve the segmentation results is validated. A agar–gel phantom has also been created, in order to test the accuracy of the parameters in conditions closer to the in vivo ultrasound imaging. This phantom contains four inclusions with different concentrations of silica. A B&K ultrasound device provides the RF data. The quantification of the segmentation quality is based on the rate of correctly classified pixels and it has been computed for all the parameters either from the field images or the phantom images. The large improvement in the segmentation results obtained reveals that the multiparametric segmentation scheme proposed in this study can be a reliable tool for the processing of noisy ultrasound data.
2005 Elsevier B.V. All rights reserved. PACS: 87.63.Df; 07.05.Pj Keywords: Ultrasonic imaging; Tissue characterisation; Image segmentation; Parametric imaging

*

Corresponding author. Tel.: +334 7243 8226; fax: +334 7243 8526. E-mail address: [email protected] (F. Davignon).

0041-624X/$ - see front matter
2005 Elsevier B.V. All rights reserved. doi:10.1016/j.ultras.2005.06.001

790

F. Davignon et al. / Ultrasonics 43 (2005) 789–801

1. Introduction Ultrasound images are of common use in medical imaging for the diagnosis of many diseases. The probe can be positioned freely by the clinician and the short acquisition time enables the possibility of getting real time data. The segmentation of these data is a very important step toward an automatic analysis and/or quantitative measurements. For example, quantitative volume parameters recovery is a unique mean of making objective reproducible and operator independent diagnosis. Thus, it is all the more important to perform a successful segmentation. Several attempts have been developed to segment ultrasound data from the echographic images: for instance, the cavities of the heart [1], liver or prostate tumours [2,3]. As ultrasound data lead to very noisy images with low contrasted regions and blurred boundaries, the conventional edge-based segmentation methods such as Canny edge detection [4] or Prewitt operators [5] are not suited to find the boundaries between tissues. For these reasons, we have focused our research on one hand on region-based segmentation methods and on the other hand on a multiparametric approach. The different tissues on an image (representing different organs or pathological and healthy tissues) are characterised by their acoustic impedances, which reflect differences in the structural organisation of the scatterers [6–9]. We have previously proposed and evaluated a multiparametric implementation of a segmentation algorithm [10,11]. The multiparametric approach consists in calculating local features from the ultrasonic data that yield to images which are maps of the feature. Then, to ensure the robustness of the method, the segmentation is not performed from the envelope image only, but from several images. This method has been completed in this study by adding acoustical parameters to the segmentation process. These parameters are representative of the scattering conditions in the tissue. They are extracted from the K-distribution (a,b) and the Nakagami distribution (m,X). These four parameters are fully described in Section 2. So, for each parameter, the computed map shows the variations of the parameter throughout the investigated area. Local measurements are thus needed for this mapping, so we must find a trade-off between the size of the local measurement window and the resolution of the parametric image. The window (either 1D or 2D) must be large enough to exhibit stationary statistics, representative of the local specificity of the tissue. Most of the parametric measurement methods implemented in the frame of this study provide qualitative measurements. Because the purpose is image segmentation, our application does not need accurate quantitative measurements but reliable relative estimations that clearly emphasise the differences between various tissues.

The implemented method for the segmentation of the multiparametric data is based on a Bayesian approach using a multiresolution Markov random field. The several parameters are computed from the image or from the radio-frequency (RF) signals. The multiparametric approach is implemented to obtain a more robust segmentation and the multiresolution approach is interesting to improve simultaneously the robustness of the segmentation and the speed of the algorithm. This paper proposes to investigate the relevance of several parameters for the segmentation of ultrasonic data. It is organised as follows. Section 2 describes the parametric images calculated from the ultrasonic data and the segmentation algorithm used. Section 3 presents on one hand the segmentation results on numerical ultrasound phantoms and, on the other hand the segmentation results on an gelatine–agar phantom. The results are discussed in Section 4 and the conclusions drawn from this study are finally summarised in Section 5.

2. Method 2.1. Acoustical parameters The conventional image used by the radiologist during an ultrasonic examination is the B-mode image computed from the envelope of the RF signal. The envelope signal contains information related to amplitude, but the frequency and phase information, related to the attenuation and to the spatial distribution of individual scatterers, is lost irretrievably by the envelope extraction operation. Several studies have proposed various methods for extracting information about the characteristics of the tissues in order to characterise the tissues and detect pathologies. Some measurements are closely related to the acoustical properties of the tissues [12], some others, like textural parameters, are not [13]. Acoustical parameters are often extracted from the RF signal and are representative of the scatterers properties within the tissues. Those parameters can provide two kinds of information. On one hand, integrated backscatter coefficients [14] and attenuation estimators [15] measure the echogenicity properties of the scatterers. On the other hand, parameters such as scatterers density, scatterers size or mean scatterer spacing (MSS) estimators characterise the spatial organisation of the scatterers [16]. Some methods propose to derive acoustical parameters from the envelope image. They are based on the numerous studies which try to link the particularities of a tissue with a specific intensity distribution on the envelope image [17]. The methods implemented in this study for the tissue characterisation are described hereafter. They lead to

F. Davignon et al. / Ultrasonics 43 (2005) 789–801

local measurements in view to provide parametric mapping that are then involved in the segmentation process. 2.2. Integrated backscatter coefficient and attenuation One of the most frequently used parameter is the integrated backscatter coefficient (IBS) [18,19]. As it is related to the acoustic impedance of the tissues, it can be useful to detect pathological cases. IBS can be computed directly from the RF signal as described in Eq. (1) [18]. X IBS ¼ PSDðf Þ ð1Þ BW

where PSD is the power spectral density, obtained from the discrete Fourier transform applied to a segment of the RF signal. The summation is performed on a limited bandwidth BW. Several approaches have been proposed in the literature to assess the ultrasonic attenuation coefficient of tissue [20,15]. The implemented method is based on the local spectral analysis and measures the central frequency shift of the RF signal [21]. As the attenuation increases with frequency, this induces a decay of the spectrum towards the low frequencies. Under the hypothesis that the power spectral density of the backscattered signal is Gaussian distributed, this decay is proportional to the attenuation and is supposed to be linear with depth. The mathematical expression of the mean central frequency is given in Eq. (2) [18]. P f  PSDðf Þ BW MCF ¼ ð2Þ IBS The estimation of these two parameters (IBS and MCF) is based on the local estimation of the power spectral density of the RF signal. This estimation is computed from 64 sampled windows. IBS and MCF are then estimated in a 20 dB bandwidth (BW) of the power spectral density. 2.3. Scatterers density The literature shows that the spatial organisation of the scatterers can be evaluated from the RF lines or from the 2D envelope image. It has been shown that these features have an influence on the intensity distribution of the envelope [22] and according to the density range and spatial organisation of the scatterers, different distribution models can be used. We propose in this study to evaluate these features from two different approaches based on the K-distribution and the Nakagami distribution. The K-distribution is a generalisation of the Rayleigh distribution. The K-distribution has two advantages over the Rayleigh distribution: it can describe the envelope statistics even when the effective

791

number of scatterers a is less than 10, and this parameter a can be estimated directly from the envelope data (Eqs. (5) and (6)). The Nakagami distribution can be considered as a generalisation of the K-distribution. It can describe the envelope statistics even when the scatterers follow a regular pattern. Its shape parameter m describes the different cases: when m is less than 1, the distribution is equivalent to the K-distribution and a can be computed from m; when m equals 1, the Nakagami distribution is equivalent to the Rayleigh distribution; when m is greater than 1, the distribution is equivalent to the rice distribution, which can describe the coherent contribution of the scatterers. This distribution can therefore describe more cases than the K-distribution. The spatial organisation of the scatterers in the tissue can be estimated by computing two parameters: a density estimator and a mean scatterer spacing estimator. The latter parameter refers to particular tissues exhibiting a regular organisation of the cells. This configuration is infrequently encountered in biological tissues and this parameter has not been considered in this study. In the case of a medium where the scatterers are not regularly spaced, the resulting echo is a coherent sum of reflected waves from randomly distributed scatterers. If the number of scatterers per resolution cell is greater than 10, the amplitude of the envelope of the signal is Rayleigh distributed [22]. For densities less than 10, the envelope is modelled with a K-distribution [22]. Its probability density function is given in Eq. (3):  a 2b b  A pðAÞ ¼ K a1 ðb  AÞ; where CðaÞ 2 rffiffiffiffiffiffiffiffiffiffiffiffi a ð3Þ b¼2 EðA2 Þ A is the envelope signal, C is the gamma real function defined so that:  Cðx þ 1Þ ¼ x  CðxÞ ð4Þ Cð1Þ ¼ 1 where a is the effective number of scatterers per resolution cell, also called effective density of scatterers. Kn is the modified second Bessel function of order n, E(A2) is the second moment of the envelope amplitude. The normalised fourth moment r4 and sixth moment r6 of the K-distribution are function of a. This leads to two estimators for the density, named a4 and a6: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 9 þ 9 þ 12r6 a6 ¼ ð5Þ a4 ¼ r6  6 r4  2 These two estimators have major drawbacks. Firstly, their relative bias and variance tend to get high when a is far from unity. Secondly, the smaller the sampling size is, the higher those two values are.

792

F. Davignon et al. / Ultrasonics 43 (2005) 789–801

So, for values far from 1 (i.e. closer to 0.1 or 10), another estimator would be better. Ossant et al. have thus proposed an estimator based on low order moments [23]. It gives good results even in the cases of a small sampling size or a density far from 1. The one used here is called a2 and its computation is based on the Eq. (6):

The right part of Eq. (6) can be computed from the image. Then, the equation can be solved numerically to obtain an estimation of a. Here, a2 was computed via a dichotomy algorithm.

the parametric images can be taken into account in the segmentation process. We assume that the observed data Y is a random field defined on a 2D rectangular grid SÆYs denotes the value of Y at the site s 2 S. A segmentation of the image into regions will be noted by X, where Xs = i means that the pixel at s belongs to region i. The number of different regions type is k. The conditional density function of Y given X, is assumed to exist and to be strictly positive and is denoted by f(y/x). The probability P(X = x) is written as P(x). The image may be segmented by estimating the pixel classification X given the observed image Y using the Maximum A posteriori estimation of X.

2.4. Nakagami distribution

^xmap ¼ arg max fP ðX ¼ x=Y ¼ yÞg

Cða þ 1Þ  CðaÞ EðA2 Þ Cð3=2Þ  ¼ 2 E ðAÞ Cð2Þ C2 ða þ 1=2Þ

ð6Þ

x

According to the scatterers density, the Nakagami distribution covers the whole range of envelope statistics [24]: Rayleigh, pre-Rayleigh (K-distribution) and postRayleigh (Rice). The probability density function (Pdf) is expressed as follows [25]: PdfðxÞ ¼

2mm x2m1 expðmx2 =XÞU ðxÞ CðmÞ  Xm

ð7Þ

where U(x) is the unit step function and C is the Gamma function. m and X can be estimated as: X ¼ EðX 2 Þ

m¼

X2 EðX 4 Þ  X2

ð8Þ

where E() is the statistical average. For the Nakagami distribution, X is the scaling parameter and m is the shape parameter. For m = 1, the Nakagami density function becomes a Rayleigh density function. The shape parameter m is closely related with the scattering conditions in the resolution cell. When the scatterers are randomly located but less than 10, then the envelope exhibits pre-Rayleigh statistics and m is lower than 1. In this case, there is a link between the Nakagami parameters and those of the K-distribution. sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2m 2m a¼ b¼2 ð9Þ 1m Xð1  mÞ When the resolution cell contains periodic scatterers, the envelope is likely to exhibits Rician statistics. In this case, m is greater than 1. So m can be a very useful parameter to estimate the scattering conditions in the resolution cell. 2.5. Segmentation algorithm The segmentation method, based on the Ashton and Parker algorithm [26] and adapted for multiparametric data [10,11] is briefly presented. We also detail how

¼ arg min f ln f ðy=xÞ  ln P ðxÞg x

ð10Þ

Hence, once the distributions of f(y/x) and P(x) are defined, the problem of segmenting an image will be reduced to a problem of minimising an energy function. We use a Markov random field to model the region process X, due to its restriction to local interaction. So, for a given neighbourhood system, the prior density P(x) can be written as a Gibbs density which has the following form [27]: ( ) X 1 P ðxÞ ¼ exp  V c ðxÞ ð11Þ Z all cliques C Here Z is a normalising constant called the partition function. Vc(x) is the clique potentials. A clique is a set of points that are neighbours of each other. In our work, we used a 6-points neighbourhood in the 3D space and consider only the two-points clique potentials which are defined as follows:   b if xs ¼ xq and s; q 2 c V c ðxÞ ¼ ; b > 0 ð12Þ þb if xs 6¼ xq and s; q 2 c The intensity of a region i at the site s is assumed to be a Gaussian distribution with mean lis and variance i ri2 s . The local class mean ls is a slowly varying function i2 of s. rs is estimated independently for each class, and is proportional to lis . Under these assumptions, lnf(y/x) may be written as: ! X

1 xs xs 2 ln f ðy=xÞ / lnðrs Þ þ y s  ls ð13Þ 2 2ðrxs s Þ s Substituting P(x) from Eq. (11) and ln f(y/x) from Eq. (13) into Eq. (10) leads to the following energy function: ! X

1 2 lnðrxs s Þ þ y s  lxs s U ðx=yÞ ¼ 2 2ðrxs s Þ s X V c ðxÞ ð14Þ þ all cliques C

F. Davignon et al. / Ultrasonics 43 (2005) 789–801

This function has two components. One constraints the region intensity to be close to the data, the other is a regularisation term which imposes a smoothness constraint. To improve the robustness of the algorithm, we consider the possibility of modifying the energy function by adding other constraints based on acoustical properties representative of each region. Let {Y1, . . . ,Yn} be a set of acoustic features calculated on each site of the volume of data. The features used must be chosen as relevant measurements of the characteristics of each region. They must be adapted to the objectives of the segmentation and to the particularities of the data processed. The parametric images are modelled in the same way than grey scale images in Eq. (13). A parametric image is assumed to be a collection of uniform or slowly varying intensities. The sharp transitions in grey level may only occur at region boundary. The parametric measurement of feature Yj at a pixel location s is denoted (yj)s. We assume that the values given by a feature Yj on a class i is x modelled by a normal distribution of mean ðmj Þs s and xs 2 variance ðrj Þ . Then the complete energy function has the following form: ! X 1 2 U ðx=y; y 1 . . . y n Þ ¼ lnðrxs s Þ þ ðy s  lxs s Þ 2ðrxs s Þ2 s þ

X

V c ðxÞ þ

all cliques C

þ

1 2  ðrxj s Þ2

ððy j Þs 

XX j

lnðrxj s Þ

s

ðmj Þxs s Þ2

! ð15Þ

Finding the global minimum of this function requires a lot of computation. Besag has proposed the Iterated Conditional Modes algorithm [28] as an alternative to simulated annealing. It modifies the label of sites in x in view to maximise the conditional density function at each site. The algorithm converges to a solution which is reasonably acceptable under a good initialisation. The initialisation is performed using the K-means clustering algorithm. K-means is a least-squares partitioning method that divide a collection of pixels into K classes. The algorithm iterates over the two following steps: (1) it computes the mean of each cluster and (2) it computes the distance of each pixel from each cluster mean and assigns each pixel to the cluster it is nearest to. The initial assignment of points to clusters is done randomly. 3. Results The relevance in the segmentation process of the various acoustic measurement performed has been evaluated through different sets of data. Texture images have first been synthesised, following either a Rayleigh distribution or a K-distribution. These

793

data were used to test the robustness and accuracy of Nakagami estimators in the ideal case where the distribution parameters are exactly known. Then, a set of simulated data has been computed using the Field software developed by Jensen [29]. Given the positions and the acoustical characteristics of the scatterers on one hand and the beam characteristics on the other hand, the software is able to compute the emitted and received wave, thus creating the corresponding RF signals. These data are used to evaluate the relevance of the calculated parameters to characterise the simulated tissues. Finally, segmentation trials have been performed on actual data acquired on a tissue mimicking phantom with a B-scanner providing the RF signals. 3.1. Distribution simulation In order to test the reliability of the estimated Nakagami parameters, we have generated texture images of size 256 · 256 pixels, following either a Rayleigh distribution, or a K-distribution. Then, the parameters m and X have been estimated on each image using 8 · 8, 16 · 16, 32 · 32 and 64 · 64 pixels estimating windows. Results showing the standard deviation of the measurements of the m parameter in the case of a Rayleigh distribution can be found in Table 1. According to the theory, the value of m should be 1 in this case. The relevance of the estimator of the Nakagami parameters applied on a synthesised K-distribution is illustrated by the Tables 2–4 in three cases corresponding to a distribution close to Rayleigh with a large number of scatterers (Table 4), close to an exponential distribution with a small number of scatterers (Table 2) and an intermediate one (Table 3). In each of these Tables, three values of the parameter b (which can be related to the size of the scatterers [30]) are presented. The results show that an accurate estimation of the Nakagami parameters is obtained on a large range of K-distribution provided that the estimation is performed on a statistically representative number of points. As the pre-Rayleigh domain corresponds to values of m between 0 and 1, we cannot afford standard deviations higher than 0.1 for reliable estimations. From Tables 1–4, we notice that windows sizes smaller than 16 · 16 pixels do not fulfil this condition. For a reliable estimation of the parameters, we choose

Table 1 Mean and standard deviation of the m parameter estimated from a simulated image following a Rayleigh distribution Window size

Mean

Standard deviation

8·8 16 · 16 32 · 32 64 · 64

1.08 1.02 1 0.99

0.25 0.13 0.06 0.03

794

F. Davignon et al. / Ultrasonics 43 (2005) 789–801

Table 2 Mean and standard deviation of the m and X parameters estimated from a simulated image following a K-distribution with a = 0.5 scatterers per resolution cell b

a

Window size

m

X

Expected value

Mean

Standard deviation

Expected value

Mean

Standard deviation

0.5

0.2

16 · 16 32 · 32 64 · 64

0.2

0.22 0.2 0.19

0.06 0.03 0.01

200

190 190 190

28 14 8

0.5

0.04

16 · 16 32 · 32 64 · 64

0.2

0.29 0.21 0.2

0.05 0.04 0.02

5000

4370 4950 4950

507 338 161

0.5

0.0133

16 · 16 32 · 32 64 · 64

0.2

0.25 0.24 0.24

0.05 0.02 0.01

45 000

44 100 44 053 44 053

5610 2530 1250

Table 3 Mean and standard deviation of the m and X parameters estimated from a simulated image following a K-distribution with a = 1 scatterer per resolution cell a

b

Window size

m Expected value

Mean

Standard deviation

1

0.2

16 · 16 32 · 32 64 · 64

0.33

0.32 0.30 0.29

0.08 0.05 0.02

200

180 180 180

19 9 3

1

0.04

16 · 16 32 · 32 64 · 64

0.33

0.35 0.33 0.32

0.08 0.04 0.02

5000

4880 4880 4880

540 240 110

1

0.0133

16 · 16 32 · 32 64 · 64

0.33

0.37 0.36 0.35

0.08 0.4 0.02

45 000

44 500 44 500 44 500

4700 2210 1100

X Expected value

Mean

Standard deviation

Table 4 Mean and standard deviation of the m and X parameters estimated from a simulated image following a K-distribution with a = 100 scatterers per resolution cell a

b

Window size

m

X

Expected value

Mean

Standard deviation

Expected value

Mean

Standard deviation

100

0.2

16 · 16 32 · 32 64 · 64

0.98

0.89 0.87 0.87

0.11 0.07 0.04

200

177 176 176

19 6 3

100

0.04

16 · 16 32 · 32 64 · 64

0.98

0.99 0.96 0.96

0.13 0.07 0.04

5000

4840 4880 4880

310 140 55

100

0.0133

16 · 16 32 · 32 64 · 64

0.98

1.00 0.98 0.97

0.13 0.06 0.03

45 000

44 600 44 600 44 600

2700 1300 600

to use estimating windows of size 32 · 32 pixels. This window size is a satisfying trade off between the tissue area considered to obtain a local measurement and an acceptable variance of the estimation. 3.2. Field simulation The numerical phantoms consist in a 60 mm depth, 30 mm wide, 10 mm thick scattering medium with

randomly disposed scatterers in order to get a density of 1.5 per resolution cell. In the Field software, the acoustic characteristics of each scatterer in the simulated tissue are defined through the parameter called ‘‘amplitude’’. This parameter controls the echogenicity of the medium. Thus, in the following of the paper, the amplitude value is referred as echogenicity. A 10 mm diameter cylindrical inclusion is placed at the centre of the phantom. The density and the echoge-

F. Davignon et al. / Ultrasonics 43 (2005) 789–801

Fig. 1. Geometrical shape of the phantom.

nicity properties are different from those of the surrounding medium. The geometry of the two-region phantom is shown in Fig. 1. The ratios between the values of density and echogenicity inside and outside the inclusion can take the following values: 0.1, 0.2, 0.5, 1, 2, 5, 10, thus giving 48 different phantoms. From this set of numerical phantoms, parameters maps have been computed using 32 · 32 pixels window size, giving IBS, density, and Nakagami parameters images. Attenuation maps based on the mean central frequency have not been computed because the field software does not take into account the frequency dependence in its attenuation model. Segmentation has been performed using the algorithm described in Section 2, with the following parameters: Number of classes: 2; Number of scale decompositions: 2; (b regularisation coefficient at the highest resolution: 0.6; Db (increment of b at each resolution): 0.2). Fig. 2 shows an example of segmentation results. The results illustrated by Fig. 2 shows clearly that the segmentation depends on the data involved in the segmentation process and using several images which represent different properties of the tissues improves the delineation of the object. When four images are used, the resulting segmentation is close to the ideal result shown in Fig. 1. The complementarity and/or the redundancy of the different data reinforce the robustness of the segmentation. An attempt to quantify the improvement observed visually has been carried out. For each phantom and each segmentation performed, we have counted the number of correctly classified pixels. The aim is to check whether adding one or more parameters improves the segmentation or not, and to what extent. Fig. 3 reports the results of this quantification. For various segmentation experiments, performed on the 48 phantoms, involving different parametric images, the percentage of correctly classified pixels is reported for different scatterers characteristics in the phantom. To begin with, some general remarks about the results presented can be made. Firstly, the different graphics exhibit an approximate U-shaped curve. This is explained by the fact that the segmentation is obviously better when the characteristics (scatterers density or echogenicity) of the two regions are highly different than when they are close. Secondly, in some cases, the segmentation with only one feature (the envelope image) shows a very large number of pixel correctly classified (larger than

795

90%). In that case, involving other parameters in the segmentation does not improve significantly the segmentation result. The cases where the density ratio is the same between the two media are put together on Fig. 3(a). We can see that the IBS image is not informative to segment regions which differ only about the echogenicity of the scatterers. Although, the echogenicity of the scatterers could be related to their size, the use of the feature b does not improve significantly the segmentation in this case. A significant effect of the feature X can be observed in the segmentation when the echogenicity ratio is larger than 0.5. The cases where the echogenicity ratio is the same between the two media are put together on Fig. 3(b). As expected, the feature a2 leads to an increase of about 10% of the correctly classified pixels for all the densities. The same feature a estimated from the Nakagami distribution parameters also improves the segmentation excepted when the density ratio is lower than 2 and it is less relevant than the feature a2. The segmentation performed with the feature X associated to the envelope is also improved when the densities of scatterers of the two media are very different. In Fig. 3(d), the two regions differ in their echogenicity (ratio equals 2). The feature a2 improves the segmentation when the density ratio is low. X is a relevant feature because it improves the segmentation when the density differences are low. Obviously, it does not bring any improvement when the segmentation is almost perfect. One can observe a very bad segmentation when the density ratio is 0.5. In this case, the density ratio is compensated by the scatterers echogenicity ratio to provides two regions which are statistically very close. The same behaviour appears on Fig. 3f, where the echogenicity ratio equals 0.5. From Fig. 3(c) and (e), where the density ratios between the two regions are respectively 2 and 0.5, we can observe that X brings the most significant improvement. The b parameter can improve the results in some cases but overall brings little improvement. From this experiment with ultrasound data simulated using the field software; we can see that for low-density tissues the segmentation can be improved by the use of adequate parameters. This study shows that the X parameter which is sensitive to the energy of the signal and the a2 parameter which measures the density of scatterers in the medium bring useful information for the segmentation. The other ones can be useful in some cases but are overall less relevant and thus should be used in very specific cases. 3.3. Gelatine–agar phantom To complete the test on simulated ultrasound data, segmentation experiments have been carried out on

796

F. Davignon et al. / Ultrasonics 43 (2005) 789–801

Fig. 2. Segmentation results of the numerical phantom with echogenicity ratio = 0.5, density ratio = 0.1. The right column shows the segmentation results involving the parametric images on the corresponding row. The first row shows a segmentation with the K-means algorithm. E: Envelope, D: Density (a2), I: IBS, X: Nakagami parameter.

ultrasonic images obtained with a home-made agar–gel phantom. The phantom is a 6 cm depth, 6 cm wide, 2.5 cm thick mass with four cylindrical inclusions. Each of these inclusions has a 1 cm diameter (see Fig. 4).

The composition of the phantom is as follows: 100 ml water, 1 g gelatine, 1 g agar–agar, 1 g of silica for the scatterers. As for the inclusions, the only difference in the composition concerns the density of scatterers. The

F. Davignon et al. / Ultrasonics 43 (2005) 789–801

(a)

(b)

(c)

(d)

(e)

(f)

797

Fig. 3. Influence of the parameters on the quality of the segmentation.

quantities of silica added in the four inclusions are 1%, 2%, 4% and 8% in volume. Acquisition of data has been made using a B&K Bscanner providing the RF signals, equipped with a 3.5 MHz probe (5 cm wide). The RF lines are converted

into digital signals by a numerical scope and stored on a computer. A software developed in the laboratory is used to control the acquisitions and to trigger them. The sampling frequency is 50 MHz and 128 RF lines have been used, thus creating a 128 · 2502 RF file.

798

F. Davignon et al. / Ultrasonics 43 (2005) 789–801

Fig. 4. Geometry of the agar–gelatine phantom.

The parameters involved in this experiment are slightly different from those used in the simulation. The a2 density parameter has not been computed because the envelope image exhibits Rayleigh statistics, thus the density estimators are irrelevant in this case. Unlike the experiments performed on simulated field data, the attenuation parameter based on the estimation of the central frequency of the RF signals can be used. The Nakagami parameters and IBS have been computed as described previously. To compute the different parameters, 14 · 342 pixels windows were used (14 RF lines and 342 samples on each line). This size corresponds to a 0.54 · 0.54 cm2 window size. The parametric images are obtained by using a sliding window with a shift between two consecutive windows of one RF line in the horizontal direction and of 24 samples in the direction of the RF line. Segmentations has been performed using the algorithm described in Section 2, with the following parameters: Number of classes: 3; Number of scale decompositions: 3; (b regularisation coefficient at the highest resolution: 0.4; Db (increment of a2 at each resolution): 0.1). Results are shown in Fig. 5. We can see that the segmentation using only the envelope image only gives very imprecise regions. The algorithm identifies the background, the four inclusions as a unique region and the shadow at the right-bottom corner of the image due to the high concentration of scatterers in the fourth inclusion. A trial with the mean central frequency parameter does not lead to any improvement, because this image exhibits very poor contrast and is mainly sensitive to the difference between the shadow and the remaining part of the image. The IBS image being quite close to the envelope image, the redundancy of the data emphasises the differences between the regions in the image, thus giving a better segmentation result. We can also observe that in this case the sizes of the regions tend to be underesti-

mated and the shapes are not very satisfying, especially for the left one which exhibits the lowest scatterers density. The X parameter is also sensitive to the energy of the backscattered signal, so the four inclusions are enhanced. Nonetheless, this feature is estimated on a rather large window, the contours of the regions are blurred and the segmentation does not separate the different areas. Combining X and IBS parameters leads to a better result: the regions are well separated and the sizes and the shapes are good, except for the left one. If we add the central frequency parameter, we can see that the differences are small: the size and the shape of the left region are better but the separation between the two regions at the left is not as good.

4. Discussion When the problem of segmenting ultrasound data is addressed, the quality of the results depends on the relevance of the data and of the segmentation tool used. Some arguments for the multiparametric approach are discussed in this section. The limitations of the implemented method are also presented. The interest of parametric imaging is to provide maps of local features which are often statistical measurements. Then, they must be calculated on windows large enough to contain statistically significant data. The experiments carried out in this study on simulated distribution show that a minimum of about 1000 independent samples (32 · 32) is required to obtain accurate and reliable measurements. The minimum area considered for local measurements on actual data depends on the resolution of the B-scanner used. In this study, the estimation windows correspond typically to 0.5 · 0.5 cm2. The relevance of the features calculated is of major importance for a good segmentation. The segmentation does not require quantitative values for the measurements performed. The relevant features are those which emphasise the relative differences between tissues. The various numerical phantoms used in our experiments have shown that the feature ‘‘scatterers density’’ is a reliable indicator for low densities (in the pre-Rayleigh domain). The features IBS and X calculated on the numerical and physical phantoms also give a good image of the scatterers density and of their echogenicity. The quality of the segmentation is evaluated through a simple method which consists in quantifying the correctly classified pixels. This method is efficient when the segmentation is evaluated on simple objects with well-known locations and shapes. In view to compare the results, all the segmentation performed either on the numerical phantom or on the physical phantom have been performed with the same

F. Davignon et al. / Ultrasonics 43 (2005) 789–801

799

Fig. 5. Segmentation results of the gelatine phantom. The right column shows the segmentation results involving the parametric images on the corresponding row. The first row shows a segmentation with the K-means algorithm. E: Envelope, I: IBS, C: mean central frequency, X: Nakagami parameter.

control parameters of the algorithm (number of classes, Markov parameters b and Db, window size). Some typical results are summarised in Table 5. It reports the percentage of correctly classified pixels obtained by using on one hand the K-means segmentation method and on the other hand the method described in

this paper applied first to the envelope image only and then with different configurations of multiparametric segmentation. The segmentation method that has been implemented, based on a markovian model of data and a multiresolution approach, is clearly more efficient that the K-means approach. The interest of the

800

F. Davignon et al. / Ultrasonics 43 (2005) 789–801

Table 5 Percentage of correctly classified pixels for different segmentation experiments performed on ultrasonic images issued from the numerical and the actual phantoms

Numerical phantom

Gel phantom

Echogenicity ratio = 0.1 Density ratio = 0.1 Echogenicity ratio = 0.1 Density ratio = 1 Echogenicity ratio = 1 Density ratio = 0.1 Echogenicity ratio = 10 Density ratio = 0.1

K-means

E

E+I

E+D

E+D+I

81.2

92.9

92.9

94.9

95

74

86.3

87.8

85.2

88.1

72.5

82

85

89.6

87.3

51.6

50.5

55.7

93.4

76.4

K-means 81.7

E 94.2

E+X 94.8

E+I+X 95.7

E+C+I+X 96.1

E: Envelope, I: IBS, D: Density (a2), X: Nakagami parameter.

multiparametric segmentation is also illustrated by an increasing rate of correctly classified pixels when parametric images are added to the envelope image in the segmentation. From the first three examples, we can see that adding parameters improves significantly the quality of the segmentation. Furthermore, we can observe that the density parameter gives better results that the IBS parameter in similar conditions. From the fourth example, we can also see that the density parameter leads to a good segmentation, whereas the IBS parameter fails to produce a significant improvement. The results on the gel phantom show that the segmentation quality is improved when several parameters are used. Involving several different data in the segmentation process, complementary to the envelope image leads to a more accurate segmentation, less sensitive to the speckle noise. Nevertheless, several limitations of the proposed methods can be pointed out. As previously said, the window size used to calculate the local information in parametric images is rather large (0.5 · 0.5 cm2 in our experiments). Consequently, the multiparametric segmentation approach does not allow to obtain better resolution on segmented image than on envelope images and then is not adapted to improve the detection of very small lesions. Although the presented results exhibit generally improvements of the segmentation using several parameters, one can remark that the choice of the parameters involved in the segmentation is of major importance and when too many parameters are considered, the quality of the results may decrease. It has also been observed on the simulated data that when the two regions exhibit close properties, the improvement of the segmentation quality is not significant in most cases. In the same way, the measurements performed on the actual data do not allow to distinguish, in the segmentation, the differences in the scatterers density of the four inclusions. The four regions are merged in the same class by the segmentation algorithm, even if the number of required classes is larger than three.

The results are also sensitive to the way parameters are computed. Window sizes and overlap used for the local measurements can modify the general shape of the parametric image and thus modify the segmentation results.

5. Conclusion The results presented in this paper show that a significant improvement of the segmentation of ultrasonic data can be achieved by the use of a multiparametric approach. This study illustrates the interest of research in tissue characterisation by ultrasound. The tissue properties extracted from the RF data or the envelope images can be either complementary or redundant information to the commonly used envelope data. The use of these parameters for the segmentation of different tissues is relatively simple because such an application does not require quantitative measurements of the tissue properties. Nevertheless, the segmentation will be significant only if the parametric images involved are relevant to emphasise the difference between the various tissues. Further work is in progress to identify the most significant parameters and to introduce a weighting of the various images in the segmentation process according to their relevance.

References [1] I. Dydenko, D. Friboulet, J.-M. Gorce, J. Dhooge, B. Bijnens, I.E. Magnin, Towards ultrasound cardiac image segmentation based on the radio frequency signal, Med. Image Anal. 7 (2003) 353–367. [2] H. Yoshida, B. Keserci, D.D. Casalino, A. Coskun, O. Ozturk, A. Savranlar, Segmentation of liver tumor in ultrasound images based on scale–space analysis of the continuous wavelet transform, presented at IEEE Ultrasonics Symposium, Sendai, Japan, 1998. [3] F.L. Lizzi, E.J. Feleppa, M. Astor, A. Kalish, Statistics of ultrasonic spectral parameters for prostate and liver examination,

F. Davignon et al. / Ultrasonics 43 (2005) 789–801

[4] [5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 44 (1997) 935– 942. J. Canny, A computational approach to edge detection, IEEE Trans. Pattern Anal. Machine Intell. 8 (1986) 679–698. J.M.S. Prewitt, Object enhancement and extraction, in: B.S.L.A. Rosenfeld (Ed.), Picture Processing and Psychopichtorics, Academic Press, New York, 1970. U. Haberkorn, I. Zuna, A. Lorentz, H. Zerban, G. Layer, G.V. Kaick, U. Rath, Echographic tissue characterisation in diffuse parenchymal liver disease: correlation of image structure with histology, Ultrason. Imaging 12 (1990) 155–170. L. Landini, L. Verrazzani, Spectral characterization of tissues microstructure by ultrasounds: a stochastic approach, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 37 (1990) 448–456. P.D. Edmonds, C.L. Mortensen, J.R. Hill, S.K. Holland, J.F. Jensen, P. Schattner, A.D. Valdes, Ultrasound tissue characterization of breast biopsy specimens, Ultrason. Imaging 13 (1991) 162–185. T.J. Hall, M.F. Insana, L.A. Harrison, G.G. Cox, Ultrasonic measurement of glomerular diameters in normal adult humans, Ultrasound Med. Biol. 22 (1996) 987. D. Boukerroui, O. Basset, A. Baskurt, G. Gimenez, A multiparametric and multiresolution segmentation algorithm of 3-D ultrasound data, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 48 (2001) 64–77. D. Boukerroui, O. Basset, A. Baskurt, A. Noble, Segmentation of ultrasound images—Multiresolution 2D and 3D algorithm based on global and local statistics, Pattern Recogn. Lett. 24 (2003) 779– 790. M.J.T.M. Cloostermans, H. Mol, J.A. Verhoef, J.M. Thijssen, In vitro estimation of acoustic parameters of the liver and correlations with histology, Ultrasound Med. Biol. 12 (1986) 39–51. F.M.J. Valckx, J.M. Thijssen, Characterization of echographic image texture by cooccurrence matrix parameters, Ultrasound Med. Biol. 23 (1997) 559–571. D. Nicholas, Evaluation of backscattering coefficients for excised tissues: results, interpretation and associated measurements, Ultrasound Med. Biol. 8 (1982) 17–28. H.A.H. Jongen, J.M. Thijssen, M.v.d. Aarssen, W.A. Verhoef, A general model for the absorption of ultrasound by biological tissues, J. Acoustic. Soc. Am. 79 (1986) 535–540. T. Varghese, K.D. Donohue, Estimating mean scatterer spacing with the frequency-smoothed spectral autocorrelation function, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 42 (1995) 451–463.

801

[17] G.E. Sleefe, P.P. Lele, Tissue characterization based on scatterer number density estimation, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 35 (1988) 749–757. [18] S. Bridal, P. Forne`s, P. Bruneval, G. Berger, Parametric (integrated backscatter and attenuation) images constructed using backscattered radio frequency signal (25–56 MHz) from human aortae in vitro, Ultrasound Med. Biol. 23 (1997) 215–229. [19] H. Rijsterborgh, F. Mastik, C.T. Lance´e, P. Verdouw, J. Roelandt, N. Bom, Ultrasound myocardial integrated backscatter signal processing: frequency domain versus time domain, Ultrasound Med. Biol. 19 (1993) 211–219. [20] R. Kuc, M. Schwartz, L.V. Micsky, Parametric estimation of the acoustic attenuation coefficient slope for soft tissue, presented at IEEE Ultrasonics Symposium, 1976. [21] T. Baldeweck, P. Laugier, A. Herment, G. Berger, Application of autoregressive spectral analysis for ultrasound attenuation estimation: interest in highly attenuating medium, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 42 (1995) 99–110. [22] B.J. Oosterveld, J.M. Thijssen, W.A. Verhoef, 3-D simulations ans experiments of the effects of diffraction and scatterer density, Ultrason. Imaging 7 (1985) 142–160. [23] F. Ossant, F. Patat, M. Lebertre, M.-L. Teriierooiterai, L. Pourcelot, Effective density estimators based on the K-distribution: interest of low and fractional order moments, Ultrason. Imaging 20 (1998) 243–259. [24] P.M. Shankar, A general statistical model for ultrasonic scattering from tissues, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 47 (2000) 727–736. [25] M. Nakagami, The m distribution—A general formula of intensity distribution in rapid fading, in: P. Press, (Ed.) Statistical Methods on Radio Wave Propagation, New York, 1960, pp. 3–36. [26] E.A. Ashton, K.J. Parker, Multiple resolution bayesian segmentation of ultrasound images, Ultrason. Imaging 17 (1995) 291– 304. [27] B.J. Besag, Spatial interaction and the statistical analysis of lattice systems, J. Royal Stat. Soc. B 26 (1974) 192–236. [28] J. Besag, On the statistical analysis of dirty pictures (with discussion), J. Royal Stat. Soc. B 48 (1986) 259–302. [29] J.A. Jensen, P. Munk, Computer phantoms for simulation ultrasound B-mode and CFL images, presented at 23rd Acoust. Imag. Symp., Boston, MA, 1997. [30] P.M. Shankar, V.A. Dumane, J.M. Reid, V. Genis, F. Forsberg, C.W. Piccoli, B.B. Goldberg, Classification of ultrasonic B-mode images of breast masses using nakagami distribution, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 48 (2001) 569–580.