International audienceUltrasound image deconvolution has been widely investigated in the literature. Among the existing approaches, the most common are based on ℓ2-norm regularization (or Tikhonov optimization) or the well-known Wiener filtering. However, the success of the Wiener filter in practical situations largely depends on the choice of the regularization hyperparameter. An appropriate choice is necessary to guarantee the balance between data fidelity and smoothness of the deconvolution result. In this paper, we revisit different approaches for automatically choosing this regularization parameter and compare them in the context of ultrasound image deconvolution via Wiener filtering. Two synthetic ultrasound images are used in order to compare the performances of the addressed methods

Basarab, Adrian

Chen, Zhouye

Kouamé, Denis

Scientific Publications of the University of Toulouse II Le Mirail

A simulation study on the choice of regularization parameter in l2-norm ultrasound image restoration

Ultrasound image deconvolution has been widely investigated in the literature. Among the existing approaches, the most common are based on ℓ2-norm regularization (or Tikhonov optimization) or the well-known Wiener filtering. However, the success of the Wiener filter in practical situations largely depends on the choice of the regularization hyperparameter. An appropriate choice is necessary to guarantee the balance between data fidelity and smoothness of the deconvolution result. In this paper, we revisit different approaches for automatically choosing this regularization parameter and compare them in the context of ultrasound image deconvolution via Wiener filtering. Two synthetic ultrasound images are used in order to compare the performances of the addressed methods

Open Archive Toulouse Archive Ouverte

     Open Archive TOULOUSE Archive Ouverte (OATAO)  OATAO is an open access repository that collects the work of Toulouse researchers and makes it freely available over the web where possible.  This is an author-deposited version published in : http://oatao.univ-toulouse.fr/ Eprints ID : 15324 The contribution was presented at EMBC 2015 :  http://embc.embs.org/2015/  To cite this version : Chen, Zhouye and Basarab, Adrian and Kouamé, Denis A simulation study on the choice of regularization parameter in l2-norm ultrasound image restoration. (2015) In: 37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC 2015), 25 August 2015 - 29 August 2015 (Milano, Italy). Any correspondence concerning this service should be sent to the repository administrator: staff-oatao@listes-diff.inp-toulouse.fr A simulation study on the choice of regularization parameter inℓ2-norm ultrasound image restoration*Zhouye Chen, Adrian Basarab, Denis Kouame´Abstract— Ultrasound image deconvolution has been widelyinvestigated in the literature. Among the existing approaches,the most common are based on ℓ2-norm regularization (orTikhonov optimization) or the well-known Wiener filtering.However, the success of the Wiener filter in practical situationslargely depends on the choice of the regularization hyperpa-rameter. An appropriate choice is necessary to guarantee thebalance between data fidelity and smoothness of the decon-volution result. In this paper, we revisit different approachesfor automatically choosing this regularization parameter andcompare them in the context of ultrasound image deconvolutionvia Wiener filtering. Two synthetic ultrasound images areused in order to compare the performances of the addressedmethods.I. INTRODUCTIONUltrasound (US) medical imaging has the advantages ofbeing noninvasive, harmless, cost-effective and portable overmany other imaging modalities such as X-ray computedtomography or Magnetic Resonance Imaging [1]. However,the limited bandwidth of the imaging transducer, the char-acteristics of the ultrasound propagation such as diffractionand the imaging system tend to degrade the quality of the USimages. To deal with this problem, significant efforts havebeen made in the last few decades.Under some weak assumptions (first order Born ap-proximation and weak scattering), the recorded US radio-frequency image can be modeled as the result of a 2Dconvolution between the tissue reflectivity function and thepoint-spread function (PSF) [2]. As a consequence, deconvo-lution methods have been intensively considered to enhancethe spatial resolution of US images. Among the existingapproaches, Wiener filtering (and variants) was one of themost explored tracks in US imaging [3], [4]. However, it iswell-known that the results of the Wiener filter or its variantslargely depends on the choice of the regularization parameter(RP) that provides a compromise between data fidelity andsmoothness of the deconvolution result. In US imaging, thechoice of the RP is either done manually (empirically or bytrial-and-error) [3], or is simply related to the signal-to-noiseratio (the RP is considered the inverse of the SNR) which isfurther estimated from the data [4]. In the general literaturerelated to Wiener filter deconvolution, several approacheshave been proposed to find an optimal value of RP in anautomatic manner. In this paper, we propose to evaluate andcompare their performances in the framework of ultrasoundimaging. To do so, two simulated US images are used.The remainder of the paper is organized as follows. Theconsidered model and the Wiener filtering approach arepresented in Section II. Section III reviews the existingapproaches of optimal RP choice. The simulation setup, thedeconvolution results and the comparison in US imaging areprovided in Section IV.II. PROBLEM FORMULATIONUnder the assumption of weak scattering and using thefirst order Born approximation, the interaction between thetissues and the propagating ultrasound waves is classicallymodeled by the following 2D convolution model [2]:y = Hx+ n (1)where y ∈ RMN×1 is the vertical concatenation of Macquired radio-frequency (RF) signals of length N (alsoknown as ”lexicographical notation” of the RF image), x ∈RMN×1 is the tissue reflectivity function using the samelexicographical notation, H ∈ RMN×MN is a block circulantwith circulant block (BCCB) matrix related to the 2D PSF ofthe system and n ∈ RMN×1 is a zero-mean additive whiteGaussian noise with variance σ2. The purpose of non-blinddeconvolution methods is to recover the tissue reflectivityfunction x from the recorded US RF image y considering thePSF known (or provided a previous estimate of the PSF). Acommon way to solve the non-blind deconvolution problemis the well-known Wiener/Tikhonov approach, that providesthe following analytical estimation of x:xˆ = (HtH + λQtQ)−1Hty = F (λ)y (2)where xˆ is the estimate of x, and Q ∈ RMN×MN is a BCCBmatrix representing the regularization operator, classicallyconsidered the identity matrix or the 2D Laplacian opera-tor. λ is a hyperparameter providing compromise betweendata fidelity and smoothness of the estimate (addressed asregularization parameter (RP) hereafter). The value of λ thatguarantees an optimal compromise between data fidelity andsmoothness depends on the variance of the noise σ2 and onthe properties of H , Q and y [5].III. METHODS FOR REGULARIZATIONPARAMETER OPTIMAL CHOICEIn the general Wiener filtering literature, several ap-proaches have been proposed to automatically fix the value ofthe RP λ [6]. We review hereafter the main approaches, thatwe further evaluate in the context of US images. Following[6], we classify them into two categories. The first includesmethods that need the knowledge of the variance of the noise(or an estimate of it) to estimate the RP, while the secondone does not have this constraint.A. Methods with knowledge of the variance of the noise1) Constrained least squares (CLS): The principle of thismethod, originally proposed in [7], is based on the residualbetween the data y and Hxˆ for a given λ, expressed as:φ(λ) =‖ y −Hxˆ ‖22=‖ y −HF (λ)y ‖22(3)The optimal λ in the sense of CLS is obtained in thecase where the residual in (3) is equal to the variance of thenoise, or more precisely to MNσ2. From this it results anequation with the unknown λ, which was shown in [7] tohave an unique solution. This optimal value of λ, that wewill denote by λCLS in this paper, is obtained by solvingthe above equation using a numerical iterative procedure.2) Degree of freedom (EDF): Similar to CLS approach,the equivalent degree of freedom method [8] estimated theRP λ from an equation relating the residual and the varianceof the noise. Based on the linear relation between y and xˆ,the residual in (3) is further expressed as:φ(λ) = σ2trace[IMN −HF (λ)] (4)where IMN ∈ RMN×MN is the identity matrix.The optimal RP given by EDF approach will be denotedhereafter by λEDF .3) Mean square error (MSE): The main idea behind thismethod is to find the value of λ that minimizes the meansquare error between x and its estimate, denoted by ǫ(λ):E(‖ ǫ(λ) ‖22) = E[‖ xˆ−x ‖22] =‖ x ‖22+E[‖ xˆ ‖22]−2E[xˆtx](5)In order to eliminate the cross-term 2E[xˆtx], which isimpossible to compute in practice because of the non-knowledge of x, it is replaced in [6] by a term depeding onthe data y and on the noise variance σ2, based on a circulantassumption in DFT domain. In this case, it is shown in [6]that minimizing (5) with espect to λ is equivalent to sovingthe following equation:‖ Q−1(IMN−HF (λ))32y ‖22= σ2trace[Q−2(IMN−HF (λ))2](6)We denote by λMSE the value of λ verifying (6) andobtained in practice with numerical optimization techniques.4) Predicitve mean square error (PMSE): Similar tothe previous approach, this method is also based on theminimization of the MSE between x and its estimate. Moreprecisely, a weighted ℓ2 error norm E[‖ Hǫ(λ) ‖22], alsocalled the predictive mean square error, is minimized [5],[9], [10]. Following the same assumptions as the previousapproach, this minimization turns out in solving the equationhereafter:E[‖ Hǫ(λ) ‖22] =‖ (IMN−HF (λ))y ‖22+2σ2trace[HF (λ)](7)The solution of (7), obtained by numerical optimization,is called λPMSE in what follows.5) Generalized Stein’s unbiased risk estimate (GSURE):The main idea behind this method is to find the RP thatminimizes the following error measure between x and itsestimate [11], [12]:e(λ) =‖ P (x− xˆ) ‖2 (8)where P = Ht(HHt)†H is the projection operator and (·)†represents the pseudo-inverse.Keeping in mind that the Wiener filter, given a value of λ,provides an analytical relation between xˆ and y, we denoteby f a function from RMN to RMN with xˆ = f(y). Asshown in [10], [11] and [12], an unbiased estimator for theprojected MSE in (8), depending on the variance of the noiseand separating the true x from its estimation f(y), is givenby:e(λ) =‖ Px ‖2 + ‖ Pf(y) ‖22−2f t(y)Ht(HHt)†y+2σ2trace(P (HtH + λQtQ)−1)(9)Minimizing (9) with respect to λ provides the optimalvalue of the RP, in the sense of GSURE approach, denotedby λGSURE .B. Methods without knowledge of the variance of the noiseThe methods presented in the previous section are allbased on the knowledge of the variance of the noise. De-pending on the application, this may be an issue and a badestimation of σ2 may cause severe errors in the optimalchoice of the RP.In the following, we briefly describe two existing ap-proaches that do not use σ2 to provide an automatical choiceof λ.1) Generalized cross-validation (GCV): Generalizedcross-validation (GCV) is one of the most popular methodsof choosing optimal RP and does not require the knowledgeof the noise variance σ2. Introduced in [5], its main idea isbased on the ”leave-one-out” principle. In the case of linearalgorithms, the GCV method leads to a simple function tominimize depending on λ:GCV (λ) =‖ (IMN −HF (λ))y ‖22[trace(IMN −HF (λ))]2(10)The value obtained by numerical minimization of (10) isdenoted by λGCV in this paper.2) Marginal likelihood (ML): This approach is based onthe minimizing of the marginal likelihood (ML) functiongiven hereafter [6]:ML(λ) =yt(IMN −HF (λ))y(det[IMN −HF (λ)])1/MN(11)Hereafter, we denote the λ minimizing (11) by λML.IV. RESULTS AND DISCUSSIONA. Experimental SettingTo evaluate the performance of these different methods, weused two synthetic ultrasound images. The two images weresimulated by 2D convolution between realistic PSFs and thetissue reflectivity functions shown as the first images in Fig.TABLE I: Deconvolution Results of Synthetic Ultrasound ImagesSNR 10dB 20dB 30dBImage Cyst Phantom Cardiac Image Cyst Phantom Cardiac Image Cyst Phantom Cardiac ImageLambda SSIM Lambda SSIM Lambda SSIM Lambda SSIM Lambda SSIM Lambda SSIM1/SNR 1,00E-01 82,05 1,00E-01 63,50 1,00E-02 86,88 1,00E-02 67,29 1,00E-03 91,61 1,00E-03 76,15CLS 2,12E-01 82,95 8,32E-01 65,32 1,60E-02 86,53 4,13E-02 67,93 1,81E-03 91,76 2,66E-03 75,24EDF 2,73E-02 76,67 3,48E-02 61,05 2,12E-03 83,74 1,91E-03 60,24 3,18E-04 89,49 2,22E-04 67,49MSE 2,37E-02 75,74 2,49E-02 40,46 2,25E-03 84,04 1,79E-03 59,69 4,09E-04 90,55 2,89E-04 70,56PMSE 3,41E-02 77,98 5,84E-02 63,73 2,02E-03 83,46 2,12E-03 61,02 2,06E-04 86,59 1,79E-04 64,27GSURE 2,37E-02 75,74 2,49E-02 40,46 2,25E-03 84,04 1,79E-03 59,69 4,09E-04 90,55 2,89E-04 70,56GCV 3,46E-02 78,06 6,05E-02 63,80 2,02E-03 83,46 2,13E-03 61,08 2,06E-04 86,59 1,79E-04 64,27ML 4,57E-02 79,46 2,79E-01 63,59 3,90E-03 86,04 1,51E-02 67,78 2,80E-04 88,65 7,00E-04 70,24Fig. 1: SNR=30dB. From left to right, the images are Cyst phantom tissue reflectivity function, Cyst phantom B-mode USimage and its deconvolution results (B-mode visualisation) for Q = I .1 and Fig. 2. Independent and identically distributed (IID)zero- mean Gaussian noise was added to the data, yieldingdifferent SNRs. This leads to the original B-mode images,as shown in Fig. 1 and Fig. 2.1) Cyst phantom image: For the first simulated image,the PSF was generated with Field II [13] and correspondsto a 3.5 MHz linear probe, sampled in the axial directionat 20 MHz. The tissue refelctivity function was obtainedby generating scatterers at uniform random positions withrandom amplitudes following gaussian distribution with zeromean and variance depending on their spatial position. Themedium consists in five hyperechoic circular cysts, fivehypoechoic circular cysts and five point reflectors.2) Simulated cardiac image: The PSF was also generatedwith Field II [13] and corresponds to a sectorial probewith the central frequency equal to 4 MHz and a axialsampling frequency of 40 MHz. The scatterer positions wereuniformly random distributed. In order to obtain an ultra-realistc simulation, the amplitudes of the scatterers wererelated to the amplitude of an in vivo cardiac image, assuggested in [14].B. Deconvolution ResultsFor both US simulations, the deconvolution results wereobtained with the Wiener filtering approach, consideringthe true PSF known. The RP was chosen with one of theapproaches in Section III or using the classical choice of λequal to the inverse of the SNR. For the methods that needthe knowledge of the SNR, the true value was employed.Two classical cases were considered using two differentregularization operators Q: the identity matrix and the 2DLaplacian operator.TableI shows the deconvolution results for each approch ofRP optimal choice. The deconvolved images are comparedto the true reflectivity function using the structural similaritymeasure (SSIM) [15]. For each simulation and for a givenSNR, we highlight in bold fonts the best result. Note that thefigures we give in Table I are obtained with a regularizationoperator equal to the identity matrix. The results follow thesame trend for a Laplacian operator.The different results may be appreciated from a quali-tative viewpoint in Fig. 1 and 2, highlighting the originalrefelctivity function, the data and the deconvolved images.All the processing was done in the RF domain. However,for visualisation reasons, we plot the corresponding B-modeimages, obtained after log-compression of the correspondingenveloppe images.C. DiscussionWe may first remark that the different strategies of RPoptimal choice provide very different values of λ and thusinfluence the quality of the deconvolution results. We shouldnote that this conclusion is less true for higher SNRs (largerthan 40 dB).Second, as predicted in [5], the results provided by PMSEand GCV approaches are very similar. Similarly, the valuesof λMSE and λGSURE are also very close. This is due tothe fact that both MSE and GSURE methods are based onthe minimization of the MSE. Despite the different ways ofminimizing the MSE, they still provide similar results.Fig. 2: SNR=30dB. From left to right, the images are Cardiac tissue reflectivity function, Cardiac B-mode US image andits deconvolution results (B-mode visualisation) for Q = Q2DL.Third, we may remark that the approaches that are notbased on the knowledge of the SNR (GCV and ML) arestill providing interesting results in US image deconvolutioncompared to the other approaches. The quantitative andqualitative results that we obtained point out that the MLapproach is more adapted to US imaging than GCV. For allthe cases, ML is among the best three methods (or very closeto the best three), with CLF and the inverse of the SNR, bothof which require the knowledge of the variance of the noise.V. CONCLUSIONSIn this paper, we have compared eight methods of choos-ing optimal RP in the context of US deconvolution viaWiener filtering. The comparison was performed on twosynthetic US images. An interesting conclusion is that theapproaches not using the knowledge of the variance of thenoise are competitive against the others and are thus veryattractive in practical situations. In a future work, theseresults should be confirmed on experimental data, combinedwith approaches of PSF and noise variance estimation.REFERENCES[1] T. L. Szabo, Diagnostic ultrasound imaging: inside out. AcademicPress, 2004.[2] J. A. Jensen, J. Mathorne, T. Gravesen, and B. Stage, “Deconvolutionof in vivo ultrasound b-mode images,” Ultrasonic Imaging, vol. 15,no. 2, pp. 122–133, 1993.[3] O. Michailovich and A. Tannenbaum, “Blind deconvolution of medicalultrasound images: A parametric inverse filtering approach,” IEEEtransactions on image processing: a publication of the IEEE SignalProcessing Society, vol. 16, no. 12, p. 3005, 2007.[4] R. Jirik and T. Taxt, “Two-dimensional blind bayesian deconvolutionof medical ultrasound images,” Ultrasonics, Ferroelectrics and Fre-quency Control, IEEE Transactions on, vol. 55, no. 10, pp. 2140–2153,2008.[5] G. H. Golub, M. Heath, and G. Wahba, “Generalized cross-validationas a method for choosing a good ridge parameter,” Technometrics,vol. 21, no. 2, pp. 215–223, 1979.[6] N. P. Galatsanos and A. K. Katsaggelos, “Methods for choosing theregularization parameter and estimating the noise variance in imagerestoration and their relation,” Image Processing, IEEE Transactionson, vol. 1, no. 3, pp. 322–336, 1992.[7] B. Hunt, “The application of constrained least squares estimation toimage restoration by digital computer,” Computers, IEEE Transactionson, vol. 100, no. 9, pp. 805–812, 1973.[8] G. Wahba, “Bayesian” confidence intervals” for the cross-validatedsmoothing spline,” Journal of the Royal Statistical Society. Series B(Methodological), pp. 133–150, 1983.[9] P. Hall and D. Titterington, “Common structure of techniques forchoosing smoothing parameters in regression problems,” Journal ofthe Royal Statistical Society. Series B (Methodological), pp. 184–198,1987.[10] S. Ramani, Z. Liu, J. Rosen, J.-F. Nielsen, and J. A. Fessler, “Regu-larization parameter selection for nonlinear iterative image restorationand mri reconstruction using gcv and sure-based methods,” ImageProcessing, IEEE Transactions on, vol. 21, no. 8, pp. 3659–3672,2012.[11] Y. C. Eldar, “Generalized sure for exponential families: Applicationsto regularization,” Signal Processing, IEEE Transactions on, vol. 57,no. 2, pp. 471–481, 2009.[12] R. Giryes, M. Elad, and Y. C. Eldar, “The projected gsure forautomatic parameter tuning in iterative shrinkage methods,” Appliedand Computational Harmonic Analysis, vol. 30, no. 3, pp. 407–422,2011.[13] J. A. Jensen, “A model for the propagation and scattering of ultrasoundin tissue,” Acoustical Society of America. Journal, vol. 89, no. 1, pp.182–190, 1991.[14] M. Alessandrini, H. Liebgott, D. Friboulet, and O. Bernard, “Simu-lation of realistic echocardiographic sequences for ground-truth vali-dation of motion estimation,” in Image Processing (ICIP), 2012 19thIEEE International Conference on. IEEE, 2012, pp. 2329–2332.[15] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Imagequality assessment: from error visibility to structural similarity,” ImageProcessing, IEEE Transactions on, vol. 13, no. 4, pp. 600–612, 2004.

http://oatao.univ-toulouse.fr/15324/1/chen_15324.pdf

A simulation study on the choice of regularization parameter in l2-norm ultrasound image restoration

Abstract

Similar works

Full text

Available Versions

Scientific Publications of the University of Toulouse II Le Mirail

Open Archive Toulouse Archive Ouverte