Bayesian inference and mathematical imaging. Part IV: mixture, - - PowerPoint PPT Presentation

Bayesian inference and mathematical imaging. Part IV: mixture, random fields, and hierarchical models.

• Dr. Marcelo Pereyra

http://www.macs.hw.ac.uk/∼mp71/

Maxwell Institute for Mathematical Sciences, Heriot-Watt University

January 2019, CIRM, Marseille.

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 0 / 55

Skin

Skin cancer is the most common form of cancer Skin melanoma kills 14000 in Europe every year Diagnosis and treatment are main public health issues

Human Skin Layers (MacNeil, 2007)

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 1 / 55

Ultrasound imaging

Ultrasound imaging US: diagnostics, routine tests, therapy and surgery US imaging of skin: new high frequency 3D ultrasound probes Study skin diseases & improve diagnosis Assess lesion boundaries prior to surgery (measure depth)

Dermis view with skin lesion outlined by the red rectangle

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 2 / 55

Ultrasound imaging

Limitations: Manual annotation of 3D images is time-consuming Strong speckle noise (SNR < 5.9dB), poor contrast & edges Segmentation is extremely operator-dependant Objective: Automatic and reliable segmentation of skin layers & lesions in 3D.

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 3 / 55

Outline

1

Statistical model for US signals (Pereyra and Batatia, 2012)

2

Supervised Bayesian US image segmentation (Pereyra et al., 2012b)

3

Unsupervised Bayesian US image segmentation (Pereyra et al., 2012a)

4

Conclusion

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 4 / 55

Outline

1

Statistical model for US signals (Pereyra and Batatia, 2012)

2

Supervised Bayesian US image segmentation (Pereyra et al., 2012b) Bayesian model Bayesian algorithm Experimental results

3

Unsupervised Bayesian US image segmentation (Pereyra et al., 2012a)

4

Conclusion

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 5 / 55

Physical Signal Model

Point scattering framework (Morse and Ingard, 1987) xn ≜ x(tn) =

M

∑

i=1

aip(tn − τi) (1) rn ≜ r(tn) = ∣

M

∑

i=1

ai [p(tn − τi) + ˜ p(t − τi)]∣ (2) M: total number of point scatterers ai: cross-section of ith scatterer τi: time of arrival of ith backscattered wave p(t) + ˜ p(t): analytic extension of the interrogating pulse p(t) Medical ultrasound Imaging: M, a and τ are unknown quantities

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 6 / 55

Statistical Signal Model

Important questions What are the possible statistical distributions of xn and rn ? What information about M, a and τ in f (xn) and f (rn) ? hola Conventional analytical answer: central limit theorem (M is very large) (Wagner et al., 1983) xn =

M

∑

i=1

aip(tn − τi) ∼ N(0,σ2

n)

rn = ∣

M

∑

i=1

ai [p(tn − τi) + ˜ p(tn − τi)]∣ ∼ Rayleigh(σn) σ2

n ∝ M⟨a2 i ⟩ is the power backscattered by the nth resolution cell

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 7 / 55

Statistical Signal Model

For many biol. tissues xn ∼ N(0,σ2) and rn ∼ Rayleigh(σ) are poor models, the empirical tails are not well modeled (Shankar et al., 1993; Shankar, 2000, 2003; Raju and Srinivasan, 2002)

Figure : Comparison of the empirical envelope pdf obtained from forearm dermis, and the corresponding estimations using the generalized gamma, Weibull and Nakagami distributions

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 8 / 55

Statistical Signal Model

How are these non-Gaussian statistics explained ? M is not large enough to enforce the CLT M is large and ai has very high variance σ2

n fluctuates strongly within homogenous regions

The signal formation model is inaccurate hola What is an appropriate non-Gaussian distribution for xn and rn ? → to answer these questions we study the limit distributions of xn and rn. Limit distribution: domain of attraction (equilibrium point) in the space of probability density functions hola The Gaussian distribution (CLT) - is only a particular case (finite variance). There are infinitely other equilibrium points.

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 9 / 55

Main results

Main results: If f (xn) converges as M → ∞ to a non-Gaussian distribution then

1 xn has a symmetric α-stable limiting distribution

xn ∼ SαS(α,γ) with α ∈ (0,2) and γ ∈ R+

2 The distribution of the scattering cross-section ai is heavy-tailed with

the same characteristic exponent α fA(ai) ∝ a−(α+1)

i

3 rn is the envelope of a complex SαS random variable

rn ∼ αRayleigh(α,γ)

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 10 / 55

Result 1: SαS statistical model

If xn converges in distribution as M → ∞, then it converges to a SαS(α,γ) distribution with α ∈ (0,2) and γ ∈ R+ hola

1 xn is a sequence of random summands aip(t − τi)

Its limit distribution must be invariant to addition

2 The characteristic function must be closed under exponentiation, only

the α-stable family has this property

3 aip(t − τi) is statistically symmetric, fX(xn) converges to a symmetric

α-stable distribution with parameters α ∈ (0,2] and γ ∈ R+

4 The case α = 2 corresponds to the Gaussian distribution and xn is

known to be not-Gaussian, we conclude that α ∈ (0,2)

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 11 / 55

Result 2: Power-law scattering cross-section distribution

Given that ai ∈ R+ and p(t − τi) ∈ [−P,P] is bounded, if xn ∼ SαS(α,γ) with α < 2, then ai follows a heavy-tailed distribution with exponent α fA(ai) ∝ a−(α+1)

i

Key idea: use necessary conditions for convergence to infer the class of fA(ai). xn is in the domain of attraction of a SαS distribution with α < 2 only if FZ(zi = aip(tn − τi)) verifies the Doebling & Gnedenko conditions (Samorodnitsky and Taqqu, 2000) C1 lim

zi→∞

FZ(−zi) 1 − FZ(zi) = C+ C− = 1 C2 lim

zi→∞

1 − FZ(zi) + FZ(−zi) 1 − FZ(lzi) + FZ(−lzi) = lα, ∀l > 0 C1 & C2 imply that FZ(zi) ∝ ∣zi∣−α + o (∣zi∣−α).

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 12 / 55

Result 2: Power-law scattering cross-section distribution

Also, zi is a product of random variables zi = aiui (Rohatgi, 1976) FZ(zi) = { ∫

−PfU(ui)[1 − FA(zi/ui)]dui

if zi < 0 FZ(0−) + ∫

P 0 fU(ui)FA(zi/ui)dui

if zi ≥ 0 with ui = p(tn − τi). Then, for zi ≫ 0 ∫

P 0 fU(ui)FA(zi/ui)dui ≈ cz−α i

This condition is verified by all power-law distributions fA(ai) = L(ai)a−(α+1)

i

where L(ai) is a slow varying function (i.e., lims→∞

L(ks) L(s) = 1)

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 13 / 55

Result 3: αRayleigh envelope distribution

The envelope rn is the amplitude of the analytic extension of xn rn cos(ϕn) = xn + yn, rn = ∣xn + yn∣ Assuming that the position of scatterers is purely random xn ∼ SαS(α,γ) ⇒ yn ∼ SαS(α,γ), yn ⊥ xn By deriving f (rn,ϕn) from f (xn,yn) and marginalizing w.r.t. ϕn rn ∼ αRayleigh(rn∣α,γ) where αRayleigh(rn∣α,γ) ≜ ∫

∞

Rayleigh(rn∣σ)S α

2 (σ2∣γ cos(πα

4 )

2 α ,1,0)dσ

= ∫

∞

rnλexp[−(γλ)α]J0(rnλ)dλ

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 14 / 55

Interpretation of α and γ

For modeling and physical interpretation purposes the scattering cross-sections can be assumed to follow a Pareto distribution fA(ai) ≅ αaα

ma−(α+1) i

am is given by am = lim

ai→∞aα i FA(ai)

hola Moreover, γ is a scale or spread parameter γ = D∗(α)

α

√ Mam where D∗(α) =

α

√

2π⟨pα

i ⟩

Γ(α) sin( πα

2 ), M is the number of scatters and ⟨pα

i ⟩ is

the α-th fractional moment of p(t − τi)

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 15 / 55

Experimental validation

Figure : Comparison of the empirical envelope pdf obtained from forearm dermis, and the corresponding estimations using the heavy-tailed Rayleigh, generalized gamma, Weibull and Nakagami distribution

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 16 / 55

Experimental validation

Comparison of distributions tails by means of a logarithmic pdfs

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 17 / 55

Outline

1

Statistical model for US signals (Pereyra and Batatia, 2012)

2

Supervised Bayesian US image segmentation (Pereyra et al., 2012b) Bayesian model Bayesian algorithm Experimental results

3

Unsupervised Bayesian US image segmentation (Pereyra et al., 2012a)

4

Conclusion

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 18 / 55

Outline

1

Statistical model for US signals (Pereyra and Batatia, 2012)

2

Supervised Bayesian US image segmentation (Pereyra et al., 2012b) Bayesian model Bayesian algorithm Experimental results

3

Unsupervised Bayesian US image segmentation (Pereyra et al., 2012a)

4

Conclusion

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 19 / 55

Segmentation Problem

r = (r1,...,rn,...,rN)T ∈ RN is a 2D or 3D B-mode image r is made up by K regions (biological tissues) C1,...,CK z = {z1,...,zN} is a label vector that maps observations r1,...,rN to tissues or classes C1,...,CK zn = k ⇔ n ∈ Ck Each region is characterized by its own α-Rayleigh statistics zn = k ⇒ rn ∼ pαR (rn∣αk,γk) We consider the maximum a posteriori (MAP) segmentation problem ˆ z = argmax

z

f (z∣r) = argmax

z

∫ ∫ f (z,α,γ∣r)dαdγ K: the number of classes (considered known), α = {α1,...,αK} and γ = {γ1,...,γK} (unknown)

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 20 / 55

Bayesian Model

We define a Bayesian model (likelihood and priors) associated to the unknown parameter vector (zT,αT,γT)T Likelihood f (r∣z,α,γ) We use an α-Rayleigh observation model f (rn∣zn = k,αk,γk) = pαR(rn∣αk,γk) Assuming observations rn are independent f (r∣z,α,γ) =

K

∏

k=1

∏

{n∣zn=k}

pαR(rn∣αk,γk) which is closely related to an α-Rayleigh mixture model f (r∣α,γ) = ∏

n K

∑

k=1

ωkfαR(rn∣αk,γk) ≈ ∑

z

f (r∣z,α,γ)

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 21 / 55

Label vector z

We consider 2D and 3D Potts Markov fields as prior distributions for z (Wu, 1982) P(z) = 1 C(β) exp ⎡ ⎢ ⎢ ⎢ ⎢ ⎣

N

∑

n=1

∑

n′∈V(n)

βδ(zn − zn′) ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ where β is the granularity coefficient, δ(⋅) is the Kronecker function and V(n) is the field’s neighborhood structure

(a) β = 0.6 (b) β = 1 (c) β = 1.2 (d) β = 1.4 Figure : Synthetic images of a 3D Potts-Markov model generated using different granularity coefficients (1 slice)

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 22 / 55

Label vector z

The proposed segmentation method uses 2D MRFs for single-slice 2D ultrasound images and 3D MRFs for multiple-slice 3D images

Figure : ultrasound image (left) and neighborhood structure V(n) (right) in 2D Figure : ultrasound image (left) and neighborhood structure V(n) (right) in 3D

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 23 / 55

α-Rayleigh parameters

characteristic index αk

Non-informative prior on αk (k = 1,...,K) αk ∼ U(0,2) this interval (0,2) covers all possible values of αk

spread γk

Inverse gamma prior on γk with hyperparameters a0 and b0 γk ∼ IG(a0,b0), k = 1,...,K where a0 = 1 and b0 = 1 to obtain a vague prior

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 24 / 55

Prior on unknown parameter vector (θ,z)

Assuming a priori independence between the parameters αk, γk and the labels z the joint prior is f (z,α,γ) = P(z)p(α)p(γ) = P(z)

K

∏

k=1

f (αk)f (γk) a0

• b0
• β
• K
• α
• γ
• z
• r

Figure : Directed acyclic graph (DAG) for the α-Rayleigh mixture model (the fixed nonrandom hyperparameters appear in dashed boxes)

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 25 / 55

Posterior density

Using Bayes theorem, the posterior distribution of (z,α,γ) is f (z,α,γ∣r) ∝ f (r∣z,α,γ)P(z)f (α)f (γ) We are interested in the marginal f (z∣r) = ∫ ∫ f (z,α,γ∣r)dαdγ MCMC Method

1 We use a Hybrid Gibbs sampler to generate samples asymptotically

distributed according to f (z,α,γ∣r)

2 We marginalize implicitly by discarding samples α(t),γ(t) 3 The samples z(t) are then used to approximate the voxel-wise MAP

estimators ˆ zn = argmaxzn f (zn∣r) for n = 1,...,N.

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 26 / 55

Outline

1

Statistical model for US signals (Pereyra and Batatia, 2012)

2

Supervised Bayesian US image segmentation (Pereyra et al., 2012b) Bayesian model Bayesian algorithm Experimental results

3

Unsupervised Bayesian US image segmentation (Pereyra et al., 2012a)

4

Conclusion

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 27 / 55

Hybrid Metropolis-within-Gibbs sampler

Generate samples asymptotically distributed according to f (z,α,γ∣r) by iteratively sampling P(z∣α,γ,r), f (α∣z,γ,r) and f (γ∣z,α,r) for t = 1,2,.. to T do — Update the label vector z — for n = 1,2,.. to N do

• 1. Draw z(t)

n

from P(zn = k∣z(t)

1∶n−1,z(t−1) n+1∶N,rn,α(t) k ,γ(t) k )

end for — Update the α-Rayleigh parameters — for k = 1,2,.. to K do

• 2. Sample α(t)

k

from f (αk∣γ(t−1),z,r)

• 3. Sample γ(t)

k

from f (γk∣α(t),z,r) end for end for Output samples z(t)

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 28 / 55

Conditional probability mass function P(z∣α,γ,r)

Generation of samples according to P[z∣α,γ,r] Sample z coordinate-by-coordinate using Gibbs moves P(zn = k∣z−n,rn,αk,γk) ∝ f (rn∣zn = k,α,γ)P(zn∣z−n) where k = 1,...,K and z−n is the vector z whose nth element has been removed. The resulting Markov random field has the following conditionals P(zn = k∣z−n,rn,αk,γk) ∝exp ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ∑

n′∈V(n)

βδ(k − zn′) ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ × rn ∫

∞

λexp[−(γkλ)αk]J0(rnλ)dλ. This step has been parallelized to reduce computing time.

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 29 / 55

Hybrid Metropolis-within-Gibbs sampler

Update α and γ

The conventional Gibbs sampler requires sampling from f (αk∣γ,z,r) and f (γk∣α,z,r). However, sampling from these distributions is not easy.

Generation using Metropolis-Hastings

Instead, we generate samples asymptotically distributed according to f (αk∣γ,z,r) and f (γk∣α,z,r) using Metropolis-Hasting moves. This results in a Metropolis-within-Gibbs sampler which also converges to the desired posterior density. Metroplis-Hastings move

• 1. Generate a sample according to an appropriate proposal distribution
• 2. Accept or reject that sample with a given probability
• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 30 / 55

Conditional probability density function p(α∣γ,z,r)

Generation of samples according to f (α∣γ,z,r) We sample α coordinate-by-coordinate using Metropolis-Hastings moves Posterior density f (αk∣γ,z,r) ∝

N

∏

{n∣zn=k}

pαR(rn∣α∗

k,γk)p(αk)

Random walk proposal α∗

k ∼ N(0,2)(α(t−1) k

,σ2

α,k)

Resulting acceptance probability 1 ∨ ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ N(0,2)(α(t−1)

k

∣α∗

k,σ2 α,k)

N(0,2)(α∗

k∣α(t−1) k

,σ2

α,k)

×

N

∏

{n∣zn=k}

pαR(rn∣α∗

k,γk)

pαR(rn∣α(t−1)

k

,γk) ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 31 / 55

Conditional probability density function p(γ∣α,z,r)

Generation of samples according to f (γ∣α,z,r) We sample γ coordinate-by-coordinate using Metropolis-Hastings moves Posterior density f (γk∣α,z,r) ∝

N

∏

{n∣zn=k}

pαR(rn∣α∗

k,γk)p(γk)

Random walk proposal γ∗

k ∼ NR+ (γ(t−1) k

,σ2

γ,k)

Resulting acceptance probability 1 ∨ ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ NR+ (γ(t−1)

k

∣γ∗

k ,σ2 γ,k)

NR+ (γ∗

k ∣γ(t−1) k

,σ2

γ,k)

×

N

∏

{n∣zn=k}

pαR(rn∣αk,γ∗

k )f (γ∗ k ∣a0,b0)

pαR(rn∣αk,γ(t−1)

k

)f (γ(t−1)

k

∣a0,b0) ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 32 / 55

Approximation of the Likelihood

Evaluating the likelihood is very time-consuming and is required at every step of the sampler and for every observation pαR(rn∣αk,γk) = rn ∫

∞

λexp[−(γkλ)αk]J0(rnλ)dλ An efficient alternative is to approximate the likelihood using the following asymptotic expansions (Sun and Han, 2008) pαR(rn∣αk,γk) =

P

∑

p=0

ap(αk,γk)r2p+1

n

+ o (r2(P+1)+1

n

) rn → 0 and pαR(rn∣αk,γk) =

P

∑

p=1

bp(αk,γk)r−αkp−1

n

+ o (r−αk(P+1)−1

n

) rn → ∞ This function has been parallelized to reduce computing time.

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 33 / 55

Outline

1

Statistical model for US signals (Pereyra and Batatia, 2012)

2

Supervised Bayesian US image segmentation (Pereyra et al., 2012b) Bayesian model Bayesian algorithm Experimental results

3

Unsupervised Bayesian US image segmentation (Pereyra et al., 2012a)

4

Conclusion

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 34 / 55

Experimental Results - Synthetic Data

Synthetic Data Generate a spatially coherent 3-component α-Rayleigh mixture with α = [1.99;1.99;1.8] and γ = [1;5;10] Estimate the posterior density f (z,α,γ∣r) using the proposed Gibbs sampler in 2D (4-pixel neighborhood, 25,000 iterations) Compute MAP estimate of z∣r and MMSE estimate of α,γ∣r

Figure : True labels, observations r (only input to the algorithm), and MAP estimates for a 3-class mixture

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 35 / 55

Experimental Results - Synthetic Data

Table : Parameter estimation

true value MMSE estimates standard deviation α1 1.99 1.99 0.002 α2 1.99 1.99 0.003 α3 1.80 1.79 0.006 γ1 1.00 1.00 0.003 γ2 5.00 5.01 0.025 γ3 10.00 9.96 0.036

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 36 / 55

Experimental Results - in vivo data

Application to real data Segmentation of an in-vivo skin lesion in a 3D high frequency ultrasound Image acquired at 100MHz with a focalized 25MHz 3D probe The number of classes K was identified visually by the clinician The granularity coefficient was fixed heuristically to β = 1 Algorithm convergence was measured using the “between-within variance criterion” (Gelman and Rubin, 1992) Results were computed using 1.000 iterations of the proposed method (900 burn-in period) Algorithm implemented in MATLAB with C-MEX and OpenMP

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 37 / 55

Experimental Results - in vivo data

(a) Dermis view with skin lesion (ROI = 450 × 200 × 16) (b) ROI (slice 7) (c) 2D Segmentation contour

Figure : Log-compressed US images of skin tumor and the estimated segmentation contours. Yellow: expert annotation, green: proposed, red: (Sarti,2005)

.

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 38 / 55

Segmentation results in 3D

(a) Slice 1 (b) Slice 3 (c) Slice 5 (d) Slice 7 (e) Slice 9 (f) Slice 11 (g) Slice 13 (h) Slice 15

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 39 / 55

Lesion reconstructed in 3D

Figure : 3D reconstruction of the melanoma tumor

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 40 / 55

Outline

1

Statistical model for US signals (Pereyra and Batatia, 2012)

2

Supervised Bayesian US image segmentation (Pereyra et al., 2012b) Bayesian model Bayesian algorithm Experimental results

3

Unsupervised Bayesian US image segmentation (Pereyra et al., 2012a)

4

Conclusion

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 41 / 55

Which β value should be used?

B-mode US ROI) MAP z (β = 0.5) MAP z (β = 0.75) MAP z (β = 1.0) MAP z (β = 1.25) MAP z (β = 1.5)

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 42 / 55

Granularity coefficient β

Previous experiments used β = 1, but β > 1 could improve results Estimate β jointly with z,α,γ from the data Inference on hierarchical Bayesian models f (z,α,γ,β∣r) Marginalize w.r.t. β: f (z∣r) = ∫ ∫ ∫ f (z,α,γ,β∣r)dαdγdβ βMAX

• a0
• b0
• β
• K
• α
• γ
• z
• r
• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 43 / 55

Granularity coefficient β

Conditional density f (β∣θ,z,r) f (β∣θ,z,r) ∝ f (r∣θ,z,β)f (θ)f (z∣β)f (β) ∝ f (z∣β)f (β) f (z∣β): Potts Markov field f (β): prior on β β ∼ U(0,βMAX) Sampling β using MH moves requires computing the ratio ratio = min{1,ξ} (3) with ξ = f (z∣β∗) f (z∣β(t−1)) f (β∗) f (β(t−1)) q(β∗∣β(t−1)) q(β(t−1)∣β∗) β∗ ∼ q(β∗∣β(t−1)) is an appropriate proposal distribution

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 44 / 55

Granularity coefficient β

Replacing f (z∣β) =

1 C(β) exp[Φβ(z)] in ξ

ξ = C(β(t−1)) C(β∗) exp[Φβ∗(z)] exp[Φβ(t−1)(z)] f (β∗) f (β(t−1)) q(β∗∣β(t−1)) q(β(t−1)∣β∗) However the ratio C(β(t−1))

C(β∗)

is intractable hola Possible solutions: Pseudo-likelihood estimators (Besag, 1975) Approximation of C(β) (Gelman and Meng, 1998; Descombes et al., 1999; Risser et al., 2009) Auxiliary variables and perfect sampling (Moller et al., 2006; Murray et al., 2006; Del Moral et al., 2006; Andrieu et al., 2010) Likelihood-free (ABC) methods (Marjoram et al., 2003; Marin et al., 2011)

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 45 / 55

Likelihood-free (ABC) Sampling

Idea: Replace f (z∣β) (intractable) by a tractable sufficient statistic η(z) f (β∣z) = f (β∣η(z))

1 Generate an auxiliary variable w ∼ PZ(w∣β) 2 Accept w if η(w) = η(z)

Indeed, η(w) = η(z) occurs with probability PZ(z∣β) hola The Gibbs potential of a Markov random fields is a sufficient statistic, i.e., η(z) =

N

∑

n=1

∑

n′∈V(n)

δ(zn − zn′) which is a scalar

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 46 / 55

Proposed Likelihood-free Metropolis Hastings Move

1: Input: {β(t−1),z(t),ν,s2

β}, number of moves M.

2: Generate β∗ ∼ N(0,B) (β(t−1),s2

β)

3: Generate w ∼ PZ(w∣β∗) through M Gibbs moves with initial state z(t) 4: if ∣η(z(t)) − η(w)∣ < ǫ then 5:

Set ratio =

f (β∗) f (β(t−1)) q(β(t−1)∣β∗) q(β∗∣β(t−1))

6:

Draw u ∼ U(0,1)

7:

if (u < ratio) then

8:

Set β(t) = β∗

9:

else

10:

Set β(t) = β(t−1)

11:

end if

12: else 13:

Set β(t) = β(t−1)

14: end if

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 47 / 55

Experimental Results - in vivo data

• Margin. β (ˆ

βMMSE = 1.02) β = 0.5 β = 0.75 β = 1.0 β = 1.25 β = 1.5

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 48 / 55

Outline

1

Statistical model for US signals (Pereyra and Batatia, 2012)

2

Supervised Bayesian US image segmentation (Pereyra et al., 2012b) Bayesian model Bayesian algorithm Experimental results

3

Unsupervised Bayesian US image segmentation (Pereyra et al., 2012a)

4

Conclusion

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 49 / 55

Conclusion

The challenges facing modern imaging sciences require a methodological paradigm shift to go beyond point estimation. In Part I we discussed how the Bayesian framework can support this paradigm shift, provided we significantly accelerate computations. In Part II we considered efficiency improvements by integrating modern stochastic and variational computation approaches. In Part III we explored methods based on convex optimisation and probability, and developed theory for MAP estimation.

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 50 / 55

In this talk we studied, though an example, Bayesian models and computation algorithms for models that are more sophisticated than the

• nes previously considered, and where deterministic approaches fail.

Thank you!

CAMM4D project Funding

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 51 / 55

Bibliography I

Andrieu, C., Doucet, A., and Holenstein, R. (2010). Particle Markov chain Monte Carlo

• methods. J. Roy. Stat. Soc. Ser. B, 72(3).

Besag, J. (1975). Statistical analysis of non-lattice data. J. Roy. Stat. Soc. Ser. D, 24(3):179–195. Del Moral, P., Doucet, A., and Jasra, A. (2006). Sequential Monte Carlo samplers. J.

• Roy. Stat. Soc. Ser. B, 68(3).

Descombes, X., Morris, R., Zerubia, J., and Berthod, M. (1999). Estimation of Markov random field prior parameters using Markov chain Monte Carlo maximum likelihood. IEEE Trans. Image Process., 8(7):945–963. Gelman, A. and Meng, X. (1998). Simulating normalizing constants: from importance sampling to bridge sampling to path sampling. Statistical Science, 13(2):163–185. Gelman, A. and Rubin, D. (1992). Inference from iterative simulation using multiple

• sequences. Stat. Sciences, 7(4):457–511.

MacNeil, S. (2007). Progress and opportunities for tissue-engineered skin. Nature, 445(1):874–880. Marin, J.-M., Pudlo, P., Robert, C. P., and Ryder, R. (2011). Approximate Bayesian Computational methods. Stat. Comput., 21(2):289–291.

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 52 / 55

Bibliography II

Marjoram, P., Molitor, J., Plagnol, V., and Tavar, S. (2003). Markov chain Monte Carlo without likelihoods. Proc. Nat. Academy Sci., 100(26):15324–15328. Moller, J., Pettitt, A. N., Reeves, R., and Berthelsen, K. K. (2006). An efficient Markov chain Monte Carlo method for distributions with intractable normalising constants. Biometrika, 93(2):451–458. Morse, P. and Ingard, K. (1987). Theoretical Acoustics. Princeton University Press, Princeton (NJ). Murray, I., Ghahramani, Z., and MacKay, D. (2006). MCMC for doubly-intractable

• distributions. In Proc. (UAI 06) 22nd Annual Conference on Uncertainty in Artificial

Intelligence, pages 359–366, Cambridge, MA, USA. Pereyra, M., Dobigeon, N., Batatia, H., and Tourneret, J.-Y. (2012a). Estimating the granularity parameter of a Potts-Markov random field within an MCMC algorithm. submitted to IEEE Trans. Image Process. Pereyra, M., Dobigeon, N., Batatia, H., and Tourneret, J.-Y. (2012b). Segmentation of skin lesions in 2D and 3D ultrasound images using a spatially coherent generalized Rayleigh mixture model. IEEE Trans. Med. Imag. to appear.

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 53 / 55

Bibliography III

Pereyra, M. A. and Batatia, H. (2012). Modeling ultrasound echoes in skin tissues using symmetric α-stable processes. IEEE Trans. Ultrason. Ferroelect. Freq. Contr., 59(1):60–72. Raju, B. I. and Srinivasan, M. A. (2002). Statistics of envelope of high-frequency ultrasonic backscatter from human skin in vivo. Ultrasonics, Ferroelectrics and Frequency Control, IEEE Transactions on, 49(7):871 –882. Risser, L., Idier, J., Ciuciu, P., and Vincent, T. (2009). Fast bilinear extrapolation of 3D Ising field partition function. Application to fMRI image analysis. In Proc. IEEE Int.

• Conf. Image Proc. (ICIP), pages 833–836, Cairo, Egypte.

Rohatgi, V. M. (1976). An Introduction to Probability Theory and Mathematical

• Statistics. Wiley-Interscience, New York (NJ).

Samorodnitsky, G. and Taqqu, M. (2000). Stable Non-Gaussian Random Processes. Chapman & Hall/CRC, New York (NJ). Shankar, P. (2000). A general statistical model for ultrasonic backscattering from

• tissues. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control,

47(3):727–736.

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 54 / 55

Bibliography IV

Shankar, P. M. (2003). A compound scattering pdf for the ultrasonic echo envelope and its relationship to k and nakagami distributions. Ultrasonics, Ferroelectrics and Frequency Control, IEEE Transactions on, 50(3):339 –343. Shankar, P. M., Reid, J. M., Ortega, H., Piccoli, C. W., and Goldberg, B. B. (1993). Use of non-rayleigh statistics for the identification of tumors in ultrasonic b-scans of the breast. Medical Imaging, IEEE Transactions on, 12(4):687 –692. Sun, Z. and Han, C. (2008). Heavy-tailed Rayleigh distribution: A new tool for the modeling of SAR amplitude images. In Proc. (IGARSS 08). IEEE Int. Geosc. and Remote Sensing Symp., volume 4, pages 1253–1256. Wagner, R. F., Smith, S. W., Sandrik, J. M., and Lopez, H. (1983). Statistics of speckle in ultrasound b-scans. Sonics and Ultrasonics, IEEE Transactions on, 30(3):156 – 163. Wu, F. Y. (1982). The Potts model. Rev. Mod. Phys., 54(1):235–268.

• M. Pereyra (MI — HWU)

Bayesian mathematical imaging 55 / 55