Preconditioned P-ULA for Joint Deconvolution-Segmentation of Ultrasound Images -- Extended Version

TL;DR

This paper introduces a novel accelerated Bayesian method using preconditioned P-ULA for joint deconvolution and segmentation of ultrasound images, significantly improving speed and quality over existing techniques.

Contribution

It proposes a new hierarchical Bayesian framework with a preconditioned P-ULA scheme and a majorization-minimization approach for faster, high-quality ultrasound image processing.

Findings

Achieves up to six times faster convergence than previous HMC methods.

Provides high-quality deconvolution and segmentation results on simulated and real ultrasound data.

Demonstrates robustness and efficiency in medical imaging applications.

Abstract

Joint deconvolution and segmentation of ultrasound images is a challenging problem in medical imaging. By adopting a hierarchical Bayesian model, we propose an accelerated Markov chain Monte Carlo scheme where the tissue reflectivity function is sampled thanks to a recently introduced proximal unadjusted Langevin algorithm. This new approach is combined with a forward-backward step and a preconditioning strategy to accelerate the convergence, and with a method based on the majorization-minimization principle to solve the inner nonconvex minimization problems. As demonstrated in numerical experiments conducted on both simulated and in vivo ultrasound images, the proposed method provides high-quality restoration and segmentation results and is up to six times faster than an existing Hamiltonian Monte Carlo method.

Tables

Table 1. TABLE I: Test settings: size of test images, data type, and availability of the ground-truth.

		Simu1	Simu2	Kidney	Thyroid	Bladder	KidneyReal
	Size	256$\times$256	256$\times$256	294$\times$354	870$\times$140	370$\times$256	350$\times$200
	Data type	Simulated	Simulated	Tissue-mimicking	Real in vivo	Real in vivo	Real in vivo
Ground-truth	TRF	✓	✓	✓	-	-	-
	$𝒢𝒢𝒟$ parameters	✓	✓	-	-	-	-
	Segmentation	✓	✓	-	-	-	-

Table 2. TABLE II: Number of iterations, computational time and MSJ per s for experiments Simu1 and Simu2.

				Iterations			Time			Mixing property
				Burn-in	Total		Duration	PP-ULA speed gain		MSJ (per s)
Simu1		P-ULA		70000	140000		2 h 27 min	12.2		665
		HMC		4000	8000		1 h 08 min	5.7		173
		PP-ULA		4000	8000		12 min	1		970
Simu2		P-ULA		70000	140000		3 h 06 min	4.8		590
		HMC		10000	20000		4 h 14 min	6.6		22
		PP-ULA		10000	20000		39 min	1		793

Table 3. TABLE III: MMSE Estimates of the noise variance and 𝒢 ​ 𝒢 ​ 𝒟 𝒢 𝒢 𝒟 \mathcal{GGD} parameters.

		Simu1						Simu2
		${\sigma }^2$	${\alpha }_1$	${\beta }_1$	${\alpha }_2$	${\beta }_2$		${\sigma }^2$	${\alpha }_1$	${\beta }_1$	${\alpha }_2$	${\beta }_2$	${\alpha }_3$	${\beta }_3$
True		0.013	1.5	1.0	0.60	1.0		33	1.5	100	1.0	50	0.50	4.0
P-ULA		0.041	2.0	0.5	0.59	1.0		122	2.0	330	2.0	3186	0.48	3.4
HMC		0.013	1.8	1.2	0.61	1.0		34	1.4	66	1.1	111	0.54	5.2
PP-ULA		0.013	1.4	0.9	0.62	1.1		35	2.3	2676	1.2	122	0.55	5.8

Table 4. TABLE IV: PSNR, SSIM, CNR and segmentation overall accuracy (OA) for simulated data.

	Simu1					Simu2
	PSNR	SSIM	CNR	OA		PSNR	SSIM	CNR	OA
Wiener - Otsu	37.1	0.57	1.26	99.5		35.4	0.63	0.97	96.0
Lasso - SLaT [29]	39.2	0.60	1.15	99.6		37.8	0.70	0.99	98.3
P-ULA	38.9	0.45	1.82	98.7		37.1	0.57	1.59	94.9
HMC	40.0	0.62	1.47	99.7		36.4	0.64	1.59	98.5
PP-ULA	40.3	0.62	1.51	99.7		38.6	0.71	1.64	98.7

Table 5. TABLE V: Number of iterations, computational time and MSJ per s for experiments on the tissue-mimicking phantom and on real data.

				Iterations			Time			Mixing property
				Burn-in	Total		Duration	PP-ULA speed gain		MSJ (per s)
Kidney		HMC		7000	14000		4 h 23 min	6.3		167
		PP-ULA		7000	14000		42 min	1		657
Thyroid		HMC		3000	6000		2 h 09 min	3.7		175
		PP-ULA		3000	6000		35 min	1		950
Bladder		HMC		5000	10000		2 h 45 min	5.2		13
		PP-ULA		5000	10000		32 min	1		1396
KidneyReal		HMC		5000	10000		1 h 49 min	5.8		11
		PP-ULA		5000	10000		19 min	1		1361

Table 6. TABLE VI: PSNR, SSIM and CNR results.

		Kidney				Thyroid		Bladder		KidneyReal
		PSNR	SSIM	CNR		CNR		CNR		CNR
Wiener		27.6	0.58	0.66		0.56		1.66		1.61
Lasso		28.5	0.59	0.67		0.99		1.76		1.76
HMC		29.5	0.62	1.10		1.52		2.23		1.88
PP-ULA		29.3	0.62	1.14		1.56		2.48		1.93

Equations

$$y=Hx+\omega ,$$

$$p(y|x,{\sigma }^2)=\frac{1}{(2\pi {\sigma }^2{)}^{n/2}}\exp \left(-\frac{\Vert y-Hx{\Vert }^2}{2{\sigma }^2}\right).$$

$$p(x|\alpha ,\beta ,z)=\prod _{i=1}^n\frac{1}{2{\beta }_{z_i}^{1/{\alpha }_{z_i}}\Gamma (1+1/{\alpha }_{z_i})}\exp \left(-\frac{|x_i{|}^{{\alpha }_{z_i}}}{{\beta }_{z_i}}\right).$$

$$p(z)=\frac{1}{C(\theta )}\exp \left({\sum }_{i=1}^n{\sum }_{j\in \mathcal{V}(i)}\theta \delta (z_i-z_j)\right),$$

$$p({\sigma }^2|y,x)\propto \mathcal{I}\mathcal{G}\left(\frac{n}{2},\frac{\Vert y-Hx{\Vert }^2}{2}\right),$$

$$p({\alpha }_k|x,\beta ,z)$$

$$p({\beta }_k|x,\alpha ,z)$$

$$p(z_i=k|x,\alpha ,\beta ,z_{\mathcal{V}(i)})\propto \frac{\exp \left({\sum }_{j\in \mathcal{V}(i)}\theta \delta (z_j-k)-\frac{|x_i{|}^{{\alpha }_k}}{{\beta }_k}\right)}{2{\beta }_k^{1/{\alpha }_k}\Gamma (1+1/{\alpha }_k)}$$

$${\mathrm{prox}}_f^Q(x)\in {\mathrm{Argmin}}_{u\in {\mathbb{R}}^n}\frac{1}{2}\Vert x-u{\Vert }_{Q^{-1}}^2+f(u).$$

$$\pi (x)=p(x|y,{\sigma }^2,\alpha ,\beta ,z)\propto \exp \left(-\frac{\Vert y-Hx{\Vert }^2}{2{\sigma }^2}-g(x)\right),$$

$${\pi }_{\gamma }(x)\propto \underset{u\in {\mathbb{R}}^n}{\sup }\pi (u)\exp \left(-\frac{\Vert u-x{\Vert }_{Q^{-1}}^2}{2\gamma }\right).$$

$$x^{(t+1)}={\mathrm{prox}}_{\gamma g}^Q({\tilde{x}}^t)+\sqrt{2\gamma }{Q}^{\frac{1}{2}}{\omega }^{(t+1)},$$

$${\tilde{x}}^{(t)}=x^{(t)}-\frac{\gamma }{{\sigma }^2}Q{H}^{\top }(H{x}^{(t)}-y).$$

$$h_{\mathcal{J}}(u,v)=\sum _{i\in {\mathcal{I}}_{\ge 1}}\frac{|u_i{|}^{{\alpha }_{z_i}}}{{\beta }_{z_i}}+\sum _{j\in \mathcal{J}}\frac{(1-{\alpha }_{z_j})v_j^{{\alpha }_{z_j}}+{\alpha }_{z_j}{v}_j^{{\alpha }_{z_j}-1}|u_j|}{{\beta }_{z_j}}.$$

$$h_{\mathcal{J}}(u,v)\ge \sum _{i\in {\mathcal{I}}_{\ge 1}\cup \mathcal{J}}\frac{|u_i{|}^{{\alpha }_{z_i}}}{{\beta }_{z_i}}=h_{\mathcal{J}}(u,(|u_i|{)}_{1\le i\le n}).$$

$$\mathrm{PSNR}=10{\log }_{10}(n{\max }_i(x_i^{\mathrm{tr}},x_i^{\mathrm{es}}{)}^2/\Vert x^{\mathrm{tr}}-x^{\mathrm{es}}{\Vert }^2).$$

$$(\forall t\in [0,+\infty [)\mathrm{d}x(t)=b(x(t))\mathrm{d}t+V(x(t))\mathrm{d}B(t),$$

$$(\forall i\in \{1,\dots ,n\})b_i(x)=\frac{1}{2}\sum _{j=1}^n{A}_{i,j}(x)\frac{\partial \log \pi (x)}{\partial x_j}+\det (A(x){)}^{\frac{1}{2}}\sum _{k=1}^n\frac{\partial }{\partial x_j}\left(A_{i,k}(x)\det (A(x){)}^{-\frac{1}{2}}\right),$$

$$(\forall t\in \mathbb{N})x^{(t+1)}=x^{(t)}+2\gamma b(x^{(t)})+\sqrt{2\gamma }{A}^{\frac{1}{2}}(x^{(t)}){\omega }^{(t)}.$$

$$x^{(t+1)}=x^{(t)}+\gamma A\nabla \log \pi (x)+\sqrt{2\gamma }{A}^{\frac{1}{2}}{\omega }^{(t+1)}.$$

$$(\forall x\in {\mathbb{R}}^n)\pi (x)=p(x|y,{\sigma }^2,\alpha ,\beta ,z)\propto \exp (-(f+g)(x)),$$

$$(\forall x\in {\mathbb{R}}^n){\pi }_{\gamma }^Q(x)=\underset{u\in {\mathbb{R}}^n}{\sup }\pi (u)\exp \left(-\frac{\Vert u-x{\Vert }_{Q^{-1}}^2}{2\gamma }\right).$$

$$x^{(t+1)}=x^{(t)}+\gamma A\nabla \log {\pi }_{\gamma }^Q(x)+\sqrt{2\gamma }{A}^{\frac{1}{2}}{\omega }^{(t+1)}.$$

$$\nabla \log {\pi }_{\gamma }^Q(x)=Q^{-1}\frac{{\mathrm{prox}}_{\gamma (f+g)}^Q(x)-x}{\gamma }.$$

$$(\forall x\in\mathbb{R}^{n})~{}~{}\log\pi^{Q}_{\gamma}(x)=-{\color[rgb]{0,0,0}(f+g)\left(\mathrm{prox}_{\gamma(f+g)}^{Q}(x)\right)-\frac{1}{2\gamma}\left\|\mathrm{prox}_{\gamma(f+g)}^{Q}(x)-x\right\|^{2}_{Q^{-1}}}.$$

$$x^{(t+1)}=x^{(t)}+\gamma A{Q}^{-1}\frac{{\mathrm{prox}}_{\gamma (f+g)}^Q(x^{(t)})-x^{(t)}}{\gamma }+\sqrt{2\gamma }{A}^{\frac{1}{2}}{\omega }^{(t+1)}.$$

$$x^{(t+1)}={\mathrm{prox}}_{\gamma (f+g)}^Q(x^{(t)})+\sqrt{2\gamma }{Q}^{\frac{1}{2}}{\omega }^{(t+1)}.$$

$$(f+g)(u)+\frac{1}{2\gamma }\Vert u-x{\Vert }_{Q^{-1}}^2=f(x)+g(u)+\frac{1}{2\gamma }\Vert u-x{\Vert }_{Q^{-1}}^2+(u-x{)}^{\top }\nabla f(x)+o(\Vert u-x\Vert ),$$

$$(f+g)(u)+\frac{1}{2\gamma }\Vert u-x{\Vert }_{Q^{-1}}^2=f(x)+g(u)+\frac{1}{2\gamma }\Vert u-x+\gamma Q\nabla f(x){\Vert }_{Q^{-1}}^2-\frac{\gamma }{2}\Vert Q^{\frac{1}{2}}\nabla f(x){\Vert }^2+o(\Vert u-x\Vert ).$$

$${\mathrm{prox}}_{\gamma (f+g)}^Q(x)$$

$$=\text{vbox}argmin{\text{displaylimits}}_{u\in {\mathbb{R}}^n}\left(g(u)+\frac{1}{2\gamma }\Vert u-x+\gamma Q\nabla f(x){\Vert }_{Q^{-1}}^2+o(\Vert u-x\Vert )\right).$$

Full text

I Introduction

In medical ultrasound (US) imaging, useful information can be drawn from the statistics of the tissue reflectivity function (TRF) to perform segmentation [1], tissue characterization [2], or classification [3]. Let $x\in\mathbb{R}^{n}$ and $y\in\mathbb{R}^{n}$ be the vectorized TRF and radio-frequency (RF) image, respectively. The following simplified model is used [4, 5]

[TABLE]

where $H\in\mathbb{R}^{n\times n}$ is a linear operator that models the convolution with the point spread function (PSF) of the probe, and $\omega\sim\mathcal{N}(0,\sigma^{2}\mathbb{I}_{n})$, with $\mathcal{N}$ the normal distribution, and $\mathbb{I}_{n}$ the identity matrix in $\mathbb{R}^{n\times n}$. This paper assumes that the PSF is known, while $\sigma^{2}>0$ is an unknown parameter to be estimated. The TRF is comprised of $K$ different tissues, which are identified by a hidden label field $z=(z_{i})_{1\leq i\leq n}\in\{1,\ldots,K\}^{n}$. For every $k\in\{1,\ldots,K\}$, the $k^{\mathrm{th}}$ region is modeled by a generalized Gaussian distribution ($\mathcal{GGD}$) [3, 6], which is parametrized by a shape parameter $\alpha_{k}\in[0,3]$, related to the scatterer concentration, and a scale parameter $\beta_{k}>0$, linked to the signal energy. Given $y$ and $H$, the aim is to estimate a deblurred image $x$ [7, 8], as well as $\sigma^{2}$, $\alpha=(\alpha_{k})_{1\leq k\leq K}$, $\beta=(\beta_{k})_{1\leq k\leq K}$, and the label field $z$. Due to the interdependence of these unknowns, it is beneficial to perform the deconvolution and segmentation tasks in a joint manner [9, 10]. This is achieved in [6] by considering a hierarchical Bayesian model, which is used within a Markov chain Monte Carlo (MCMC) method [11] to sample $x$, $\sigma^{2}$, $\alpha$, $\beta$, and $z$ according to the full conditional distributions. Despite promising results in image restoration and segmentation, the method in [6] is of significant computational complexity, in particular due to the adjusted Hamiltonian Monte Carlo (HMC) method [12, 13] used to sample the TRF. Recently, efficient and reliable stochastic sampling strategies have been devised [14, 15, 16] using the proximity operator [17], which is known as a useful tool for large-scale nonsmooth optimization [18]. In this work, we investigate an MCMC algorithm to perform the joint deconvolution and segmentation of US images, where the TRF is sampled with a scheme inspired from the proximal unadjusted Langevin algorithm (P-ULA) [15]. The latter generates samples according to an approximation of the target distribution without acceptance test, while being geometrically ergodic whereas classical unadjusted Langevin algorithms may have convergence issues.

I-A Main contributions

Our contributions include i) the proposition of an original accelerated preconditioned version of P-ULA (PP-ULA), which relies on the use of a variable metric forward-backward strategy [19, 20], ii) an efficient solver based on the majorization-minimization (MM) principle to tackle the involved nonconvex priors, and iii) a new hybrid Gibbs sampler yielding a substantial reduction of the computational time needed to perform joint high-quality deconvolution and segmentation of both simulated and in vivo US images.

This article is organized as follows: Section II describes the investigated Bayesian model and sampling strategy. Section III focuses on the proposed TRF sampling method. Numerical experiments are finally presented in Section IV.

II Bayesian Model

II-A Priors

Fig. 1 illustrates the hierarchical model used to perform a joint deconvolution-segmentation of ultrasound images. The following likelihood function is derived from (1)

[TABLE]

The TRF is a mixture of $\mathcal{GGD}$s which, under the assumption that the pixel values are independent given $z$, leads to

[TABLE]

Uninformative Jeffreys priors are assigned to the noise variance and scale parameters, while the shape parameters are assumed to be uniformly distributed between 0 and 3. The labels $z$ are modeled by a Potts Markov random field with prior

[TABLE]

with $\delta$ the Kronecker function, $C(\theta)>0$ a normalizing constant, $\theta>0$ a granularity coefficient, and $\mathcal{V}(i)$ the set of four closest neighbours of the $i^{\mathrm{th}}$ pixel.

II-B Conditional distributions

The different variables are sampled according to their conditional distributions, which are provided in this section. The conditional distribution of the noise variance is derived from the Bayes theorem as follows

[TABLE]

where $\mathcal{IG}$ denotes the inverse gamma distribution. Assuming that the different regions have independent shape and scale parameters for every $k\in\{1,\ldots,K\}$, we obtain

[TABLE]

with $\mathcal{I}_{k}=\{i\in\{1,\ldots,n\}\mid z_{i}=k\}$, $n_{k}$ the number of elements in $\mathcal{I}_{k}$, and $\mathbf{1}_{[0,3]}$ the characteristic function of $[0,3]$. Samples for $\alpha_{k}$ are drawn from (6) by using a Metropolis-Hastings (MH) random walk. For every pixel $i\in\{1,\ldots,n\}$ and every region $k\in\{1,\ldots,K\}$, the Bayes rule applied to the segmentation labels leads to

[TABLE]

where $z_{\mathcal{V}(i)}$ denotes the label values in the neighborhood of $i$. As a consequence, the label $z_{i}$ is drawn from $\{1,\ldots,K\}$ using the above probabilities (suitably normalized).

III Preconditioned P-ULA

III-A Notation

Let $\mathcal{I}_{<1}=\{i\in\{1,\ldots,n\}\mid\alpha_{z_{i}}<1\}$ and $\mathcal{I}_{\geq 1}=\{1,\ldots,n\}\setminus\mathcal{I}_{<1}$. Let $\mathcal{S}_{n}$ denote the set of symmetric positive definite matrices in $\mathbb{R}^{n\times n}$, and let ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ denote the spectral norm. For every $Q\in\mathcal{S}_{n}$, let $\|\cdot\|_{Q}=\left\langle\cdot,Q\cdot\right\rangle^{1/2}$. For every function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\{+\infty\}$, the proximity operator of $f$ at $x\in\mathbb{R}^{n}$ with respect to the norm induced by $Q^{-1}\in\mathcal{S}_{n}$ is defined as follows [17],

[TABLE]

If $Q$ is not specified, then $Q=\mathbb{I}_{n}$. If $\mathrm{prox}_{f}$ is simple to compute, then the solution to (9) for an arbitrary $Q\in\mathcal{S}_{n}$ can be obtained by using the dual forward-backward (DFB) algorithm [21], summarized in Algorithm 1. If $f$ is proper, lower semicontinuous, and convex, then the sequence $(u^{(p)})_{p\in\mathbb{N}}$ generated by Algorithm 1 converges to $\mathrm{prox}_{f}^{Q}(x)$.

III-B Sampling the TRF

The conditional distribution of the TRF is

[TABLE]

where $(\forall x\in\mathbb{R}^{n})$ $g(x)=\sum_{i=1}^{n}\beta_{z_{i}}^{-1}|x_{i}|^{\alpha_{z_{i}}}$. Let $\gamma>0$ and let $Q\in\mathcal{S}_{n}$ be a preconditioning matrix used to accelerate the sampler [22]. Following [15], $\pi(x)$ is approximated by

[TABLE]

As shown in Appendix, the Euler discretization of the Langevin diffusion equation [23] applied to $\pi_{\gamma}$ with stepsize $2\gamma$ and preconditioning matrix $Q$ leads to

[TABLE]

where $\omega^{(t+1)}\sim\mathcal{N}(0,\mathbb{I}_{n})$ and

[TABLE]

Since the proposed sampling strategy is unadjusted, (12) is not followed by an acceptance test. The bias with respect to $\pi$ increases with $\gamma$, as the speed of convergence of the algorithm. A compromise must be found when setting $\gamma$.

When $\mathcal{I}_{<1}$ is not empty, we use the MM principle [24] to replace the nonconvex minimization problem involved in the computation of $\mathrm{prox}_{\gamma g}^{Q}$ with a sequence of convex surrogate problems. Let $\mathcal{J}\subset\mathcal{I}_{<1}$. We define $h_{\mathcal{J}}$ at every $(u,v)\in\mathbb{R}^{n}\times\mathbb{R}_{+*}^{n}$ by

[TABLE]

From concavity, we deduce that, for every $v\in\mathbb{R}_{+*}^{n}$ and $u\in\mathbb{R}^{n}$ such that $\mathcal{J}\subset\{i\in\mathcal{I}_{<1}\mid|u_{i}|>0\}$, the following majoration property holds

[TABLE]

Since $h_{\mathcal{J}}(\cdot,v)$ is convex and separable, its proximity operator in the Euclidean metric is straightforward to compute. More precisely, for every $i\in\mathcal{I}_{\geq 1}$, $\eta>0$ and $s\in\mathbb{R}$, $\mathrm{prox}_{\eta^{-1}|\cdot|^{\alpha_{z_{i}}}}(s)$ has either a closed form [25] or can be found using a bisection search in $[0,|s|]$. Algorithm 1 can then be called, in order to compute the proximity operator of $h_{\mathcal{J}}(\cdot,v)$ in any metric $Q\in\mathcal{S}_{n}$. This leads to Algorithm 2 which generates a sequence $(u^{(q)})_{q\in\mathbb{N}}$ estimating $\mathrm{prox}_{\gamma g}^{Q}(\tilde{x}^{(t)})$.

The resulting Gibbs sampler is summarized in Algorithm 3.

IV Numerical experiments

IV-A Experimental settings

Six experiments are presented. Simu1 and Simu2 refer to simulated images with two and three regions, respectively. Kidney denotes the tissue-mimicking phantom produced from $10^{6}$ scatterers uniformly distributed over a digital image of human kidney tissue provided with the Field II ultrasound simulator [26]. The amplitude of each scatterer is produced using a zero-mean Gaussian distribution whose variance is linked to the amplitude of the point on the digital image. The PSF for the aforementioned simulations is obtained with Field II and corresponds to a 3.5 MHz linear probe. We also perform tests on three real ultrasound images. Thyroid denotes a real RF image of thyroidal flux obtained in vivo with a 7.8 MHz probe. The unknown PSF is identified using the RF image of a wire cross-section which was acquired with the same probe. Since the diameter of the wire is of the order of a few $\mu$m, its cross-section can almost be viewed as a point. Thus, its RF image provides a good approximation of the PSF. Finally, Bladder and KidneyReal refer to the RF images of a mouse bladder and mouse kidney, respectively. Both images were obtained in vivo with a 20 MHz probe. The PSF for these two real images is estimated using the same method as for Thyroid. The number of regions $K$ is set to 2 for Simu1 and KidneyReal, and it is set to 3 for Simu2, Kidney, Thyroid and Bladder. The test settings and images can be found in Table I and Figs. 4 and 6 (first column), respectively.

The TRF is initialized using a pre-deconvolved image obtained with a Wiener filter, while the segmentation is initialized by applying a $7\times 7$ median filter and the Otsu method [27] to the B-mode of the initial TRF. Shape and scale parameters are randomly selected in $[0.5,1.5]$, and $[1,200]$, respectively. The granularity parameter $\theta$ for the Potts model (4) is adjusted to ensure that the percentage of isolated points in the segmentation, obtained with a $3\times 3$ median filter, is close to 0.05, 0.1, 0.8, 0.08, 0.08 and 0.08 for Simu1, Simu2, Kidney, Thyroid, Bladder and KidneyReal, respectively.

IV-B Comparisons and evaluation metrics

All computational times are given for simulations run on Matlab 2018b on an Intel Xeon CPU E5-1650 3.20GHz. The code for the proposed method is available online111https://github.com/mccorbineau/PP-ULA. In addition to comparing Algorithm 3 with HMC [6], the quality of the deconvolution is compared with the one obtained with a Wiener filter, where the noise level has been estimated as in [28], and with the solution to the Lasso problem, where the regularization weight is set i) manually when the ground-truth is not available, or ii) using a golden-section search to maximize the peak signal-to-noise ratio (PSNR) defined as (with $x^{\rm tr}$ the true TRF and $x^{\rm es}$ the estimated one)

[TABLE]

We also compare our results with the segmentation given by Otsu’s method [27] applied to the Wiener-deconvolved image, and with the SLaT method [29] applied to the Lasso-deconvolved image. PP-ULA is used with $\gamma=0.09$ and $Q$ an approximation of the inverse of the Hessian of the differentiable term in (10) [30], $Q=\sigma^{2}(H^{\top}H+\lambda\mathbb{I}_{n})^{-1}$, with $\lambda=0.1$ so that $Q$ is well-defined. We have also computed the structural similarity measure (SSIM) [31] of the restored TRF and the contrast-to-noise ratio (CNR) [32] between two windows from different regions of the B-mode TRF images. The segmentation is evaluated according to the percentage of correctly predicted labels, or overall accuracy (OA). The minimum mean square error (MMSE) estimators of all parameters in HMC and PP-ULA are computed after the burn-in regime. Moreover, to evaluate the mixing property of the Markov chain after convergence, we compute the mean square jump (MSJ) per second, which is the ratio of the MSJ to the time per iteration. The MSJ is obtained using $T$ samples of the TRF $(x^{t_{0}+1},\ldots,x^{t_{0}+T})$ generated after the burn-in period, i.e. $\mathrm{MSJ}=\left(\frac{1}{T-1}\sum\nolimits_{t=1}^{T-1}\|x^{(t_{0}+t)}-x^{(t_{0}+t+1)}\|^{2}\right)^{1/2}.$

IV-C Results on simulated data

The convergence speed of Algorithm 3 is empirically observed for Simu1 and Simu2, as illustrated in Fig. 2, where we also display the results of the non-preconditioned P-ULA, for which $Q=\mathbb{I}_{n}$ and $\gamma=1.99\sigma^{2}/{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|H\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}$. Comparing P-ULA and PP-ULA on these simulated data allows us to study the effect of adding a preconditioner in the proposed sampling scheme. As reported in Table II, P-ULA needs more iterations and more time to converge than PP-ULA: the proposed method is 12.2 and 4.8 times faster than P-ULA on Simu1 and Simu2, respectively. In addition, from Table III and Fig. 3, we deduce that P-ULA is more biased than PP-ULA, which samples correctly the target distributions. Finally, as one can see in Fig. 2 and Table IV, P-ULA leads to lower PSNR, SSIM and OA values than PP-ULA. These results clearly emphasize the benefits of preconditioning in this example.

From Table II, PP-ULA is 5.7 and 6.6 times faster than HMC on Simu1 and Simu2 and has better mixing properties, as shown by the MSJ per second. Visual results from Fig. 4 and CNR values in Table IV show that the contrast obtained with PP-ULA is better than with competitors on Simu2, and is second best after P-ULA on Simu1. However, it should be noted that the PSNR and SSIM obtained on Simu1 with P-ULA are much lower than with the other methods. In addition, the PSNR and SSIM values from Table IV obtained with PP-ULA are equivalent or higher than all competitors for these two experiments. Visual segmentation results are shown in Fig. 5, and OA values can be found in Table IV. For these simulated images, more pixels are correctly labeled with PP-ULA than with competitors.

IV-D Results on a tissue-mimicking phantom and on real data

The convergence of Algorithm 3 is also empirically observed for the experiments on the tissue-mimicking phantom and on real data, i.e. Kidney, Thyroid, Bladder and KidneyReal. As mentioned in Table V, the proposed method leads to a significant acceleration since it is between 3.7 and 6.3 times faster than HMC on these experiments. Visual results from Fig. 6 and CNR values in Table VI show that the contrast obtained with PP-ULA is better than with competitors on all these test images. In addition, the PSNR and SSIM values from Table VI obtained with PP-ULA on the Kidney experiment are equivalent or higher than all competitors. Although the ground-truth of the segmentation is not available for these experiments, one can see from the visual segmentation results shown in Fig. 7, that the segmentation based on the Potts model (PP-ULA and HMC) gives more homogeneous areas than Otsu, and recovers more details than SLaT.

V Conclusion

We investigated a new method based on a preconditioned proximal unadjusted Langevin algorithm for the joint restoration and segmentation of ultrasound images, which showed faster convergence than an existing Hamiltonian Monte Carlo algorithm. Hence, the proposed method has the potential to speed-up the approach proposed in [1] for the segmentation of ultrasound images. Another direction for future work is to extend this framework to a spatially variant, possibly unknown, PSF.

Appendix

In this section, after reminding results about the Langevin diffusion and its discretization using Euler’s scheme, we provide details about the derivation of the proposed method (12) used to sample the TRF.

V-A Discrete Langevin diffusion

An $n$-dimensional Langevin diffusion is a continuous time Markov process $(x(t))_{t\in[0,+\infty[}$ taking its values in $\mathbb{R}^{n}$, which is the solution to the following stochastic differential equation [23],

[TABLE]

where $(B(t))_{t\in[0,+\infty[}$ is a Brownian motion with values in $\mathbb{R}^{n}$, and for every $x\in\mathbb{R}^{n}$, $V(x)\in\mathbb{R}^{n\times n}$ is the volatility matrix and $b(x)=(b_{i}(x))_{1\leq i\leq n}\in\mathbb{R}^{n}$ is the drift term defined as

[TABLE]

where $A(x)=V(x)V(x)^{\top}=(A_{i,j}(x))_{1\leq i,j\leq n}$ is a symmetric positive definite matrix, $\mathrm{det}(A(x))$ denotes its determinant, and $\pi$ is the density of the stationary distribution of the diffusion. Here, we take $(\forall x\in\mathbb{R}^{n})$ $\pi(x)=p(x|y,\sigma^{2},\alpha,\beta,z)$ defined in (10). Euler’s discretization scheme applied to (15) leads to the following target posterior distribution, which can be used to generate a Langevin Markov chain.

[TABLE]

Hereabove, $\omega^{(t)}\sim\mathcal{N}(0,\mathbb{I}_{n})$ and $\gamma>0$ is the discretization stepsize that controls the length of the jumps, while the scale matrix $A(\cdot)$ drives their direction. Instead of taking $A(\cdot)=\mathbb{I}_{n}$ as in the standard Metropolis adjusted Langevin algorithm, we follow [22, 19] and use a preconditioning matrix $A$ to accelerate the Langevin scheme, which leads to

[TABLE]

V-B Approximation of the target diffusion

For every $x\in\mathbb{R}^{n}$, let $f(x)=\|y-Hx\|^{2}/(2\sigma^{2})$. From (10), the target distribution statisfies the following relation,

[TABLE]

where $(\forall x\in\mathbb{R}^{n})$ $g(x)=\sum_{i=1}^{n}\beta_{z_{i}}^{-1}|x_{i}|^{\alpha_{z_{i}}}$. Let $\gamma>0$ and $Q\in\mathcal{S}_{n}$. Following [15], we replace $\pi$ by its Moreau approximation $\pi_{\gamma}^{Q}$ defined in (11) and recalled below,

[TABLE]

Note that we dropped the normalization constant. Moreover, at the difference of [15], we introduce the preconditioning matrix $Q$ for convergence acceleration purposes. When $Q$ is not specified, the identity matrix is used, i.e. $Q=\mathbb{I}_{n}$. Hence, the approximated version of (18) reads

[TABLE]

We can then deduce the following result when $g$ is convex.

Proposition 1

For every $\gamma>0$, $Q\in\mathcal{S}_{n}$ and $x\in\mathbb{R}^{n}$, if $(\forall k\in\{1,\ldots,K\})$ $\alpha_{k}\geq 1$, then we have

[TABLE]

Proof:

By definition of $\pi^{Q}_{\gamma}$, we have

[TABLE]

Hence, applying [33, Lemma 2.5] in the metric induced by $Q^{-1}$ directly leads to the result. ∎

From Proposition 22, (21) becomes

[TABLE]

It can be noted that, in Proposition 22, $g$ is assumed to be convex, which is not necessarily satified in our case. However, for simplicity, we take the discrete scheme (24) even in the nonconvex case. Finally, we take $A=Q$, which leads to

[TABLE]

V-C Forward-backward approximation

By definition, $f$ is differentiable on $\mathbb{R}^{n}$ and its gradient $\nabla f=H^{\top}(H\cdot-y)/\sigma^{2}$ is Lipschitz-continuous on $\mathbb{R}^{n}$. It is worth noting that the computation of the proximity operator of the sum of two functions is generally intractable [34]. Hence, as suggested in [15], we use a first-order Taylor expansion to approximate the proximity operator of $f+g$ and introduce a forward step in PP-ULA iteration. Let $o$ denotes Landau’s notation.222Following Landau’s notation, we will write that $F(u)=o(\|u-x\|)$, where $F:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and $x\in\mathbb{R}^{n}$, if $F(u)/\|u-x\|\rightarrow 0$ as $u\rightarrow x$. Let $x\in\mathbb{R}^{n}$, using $(\forall u\in\mathbb{R}^{n})$ $f(u)=f(x)+(u-x)^{\top}\nabla f(x)+o(\|u-x\|)$, we have

[TABLE]

which can be re-written as

[TABLE]

Hence, the proximity operator of $f+g$ can be expressed as follows,

[TABLE]

In addition, we have

[TABLE]

Therefore, when $\gamma$ is small, $\mathrm{prox}_{\gamma g}^{Q}(x-\gamma Q\nabla f(x))$ is a good approximation of $\mathrm{prox}_{\gamma(f+g)}^{Q}(x)$. Plugging this preconditioned forward-backward scheme [18] in (25) leads to the proposed sampling method

[TABLE]