Total variation denoising with iterated conditional expectation - - PowerPoint PPT Presentation

TV-ICE denoising Other ICE tasks Discussion

Total variation denoising with iterated conditional expectation

C´ ecile Louchet1 and Lionel Moisan2

1 Universit´

e d’Orl´ eans, Institut Denis Poisson, France

2 Universit´

e Paris Descartes, MAP5, France

Workshop IHP “Statistical modeling for shapes and imaging” March, 12th 2019

TV-ICE denoising Other ICE tasks Discussion

TV restoration of images

Image formation model v = Au + n v ∈ RΩ′: observed image A : RΩ → RΩ′: linear operator (A = Id → denoising; A = k ∗ · → deblurring...) n: Gaussian additive white noise ∼ N(0, σ2) u ∈ RΩ: image that we want to estimate. Rudin-Osher-Fatemi image recovery Choose ˆ uROF = arg min

u∈RΩ E(u) := Au − v2 + λTV (u)

Total Variation: TV (u) = ∇u1 λ ≥ 0 is a user-controlled regularity parameter.

TV-ICE denoising Other ICE tasks Discussion

TV restoration of images

Bayesian viewpoint ˆ uROF is a Maximum A Posteriori in a Bayes framework: ˆ uROF = arg min

u Au − v2 + λTV (u)

= arg max

u

1 Z e− Au−v2

2σ2

e−βTV (u) (where β = λ 2σ2 ) = arg max

u

P(v|u) P(u) = arg max

u

P(u|v) with    P(v|u) = 1

Z e− Au−v2

2σ2

(image formation model) P(u) =

1 Z ′ e−βTV (u)

(prior distribution)

TV-ICE denoising Other ICE tasks Discussion

Restoration with TV-LSE

We have ˆ uROF = arg maxu P(u|v). Definition: image restoration by TV-Least Square Estimator [1] ˆ uLSE = E[u|v] = 1 Z

• RΩ u e−

1 2σ2 (Au−v2+λTV (u)) du

No staircasing in LSE denoising (A = Id) ∀x, y ∈ Ω, the set {v ∈ RΩ : ˆ uLSE(x) = ˆ uLSE(y)} has measure 0. Computation of TV-LSE For each x ∈ Ω, ˆ uLSE(x) = 1 Z

• RΩ u(x) e−

1 2σ2 E(u) du.

integral on RΩ where |Ω| = number of pixels ≈ 106... MCMC techniques but with convergence in O(1/ √ N).

[1] L., Moisan (2013). Posterior expectation of the total variation model: Properties and experiments. SIIMS.

TV-ICE denoising Other ICE tasks Discussion

Outline

1

TV denoising with Iterated Conditional Expectations

2

Other (imaging?) tasks with ICE Deblurring and inverse problems regularized with TV TV-ICE denoising for Poisson noise ICE of a convex functional ICE of a convex set

TV-ICE denoising Other ICE tasks Discussion

Outline

1

TV denoising with Iterated Conditional Expectations

2

Other (imaging?) tasks with ICE Deblurring and inverse problems regularized with TV TV-ICE denoising for Poisson noise ICE of a convex functional ICE of a convex set

TV-ICE denoising Other ICE tasks Discussion

The idea of TV-ICE denoising

Recall in the case A = Id: ˆ uLSE(x) =

• RΩ u(x) e− u−v2+λTV (u)

2σ2

du

• RΩ e− u−v2+λTV (u)

2σ2

du Idea: integrating one variable at a time

• R u(x) e− u−v2+λTV (u)

2σ2

du(x)

• R e− u−v2+λTV (u)

2σ2

du(x) =: fv(x)(u(Nx)) This is the posterior expectation of u(x) conditionally to u(xc). It depends on the values of u(xc). But we can iterate: convergence hopefully? From now on: Nx = 4-neighbor system and TV (u) = 1 2

• x∈Ω
• y∈Nx

|u(y) − u(x)|.

TV-ICE denoising Other ICE tasks Discussion

One iteration: closed formula

If u(Nx) = {a, b, c, d} with a ≤ b ≤ c ≤ d and if v(x) = t, then For any n ∈ N∗, for any sorted n-uple (aj) and for any n-uple (βj), we have ft(u(Nx)) = t − n

i=0 µiI t µi,νi(ai, ai+1)

n

i=0 I t µi,νi(ai, ai+1)

where {a1, . . . , a4} = {u(Nx)} and −∞ = a0 ≤ a1 ≤ · · · ≤ a5 = +∞, ∀i, µi = λ 2

n

• j=1

εi,j, νi = −λ

n

• j=1

εi,jaj εi,j =

• 1

if i ≥ j −1

• therwise

and I t

µ,ν(a, b) =

• erf

b − t + µ σ √ 2

• − erf

a − t + µ σ √ 2

• e−

1 2σ2 (2µt−µ2+ν).

TV-ICE denoising Other ICE tasks Discussion

Theorem and definition of TV-ICE Consider an image v : Ω → R and λ, σ > 0. The sequence (un)n≥0 defined recursively by u0 and ∀x ∈ Ω, un+1(x) = fv(x)(un(Nx)) converges linearly to an image ˆ uICE independent of u0 and satisfies ∀x ∈ Ω, ˆ uICE(x) = Eu|v[u(x) | u(xc) = ˆ uICE(xc)]. Idea of the proof: we define Fv by un+1 = Fv(un). Then Fv(u)(x) = fv(x)(u(Nx)). Fv is C1 and monotone: w1 ≤ w2 ⇒ Fv(w1) ≤ Fv(w2) ft−c(w(Nx) − c) = ft(w(Nx)) − c and implies Jac Fv∞ < 1 Kw =

• min(minΩ v, minΩ w), max(maxΩ v, maxΩ w)

Ω satisfies Fv(Kw) ⊂ Kw for any w.

TV-ICE denoising Other ICE tasks Discussion

Properties of TV-ICE denoising

ICE is not LSE. Proof: LSE is a prox, ICE is not. No staircasing Let x and y be neighbor pixels. The set {v ∈ RΩ : ˆ uICE(x) = ˆ uICE(y)} has measure 0. Proof: v → ˆ uICE is C1. Recovery of edges v(x) − 2λ ≤ ˆ uICE(x) ≤ v(x) + 2λ. → a strong local contrast essentially persists. Proof: ft(a, b, c, d) is a weighted average (with positive coefficients) of t, t ± λ, and t ± 2λ, it belongs to [t − 2λ, t + 2λ]. This latter property is shared with ˆ uROF and ˆ uLSE.

TV-ICE denoising Other ICE tasks Discussion

noisy ROF ICE LSE

TV-ICE denoising Other ICE tasks Discussion

noisy ROF ICE LSE

TV-ICE denoising Other ICE tasks Discussion

noisy ROF ICE LSE

TV-ICE denoising Other ICE tasks Discussion

noisy ROF ICE LSE

TV-ICE denoising Other ICE tasks Discussion

Convergence curves for different algorithms of TV-denoising

1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1 100 100 200 300 400 500 600 700 800 900 1000 TV-LSE ROF-dual ROF-primal-dual ROF-Weiss-Nesterov Huber-ROF TV-ICE

TV-ICE denoising Other ICE tasks Discussion

Outline

1

TV denoising with Iterated Conditional Expectations

2

Other (imaging?) tasks with ICE Deblurring and inverse problems regularized with TV TV-ICE denoising for Poisson noise ICE of a convex functional ICE of a convex set

TV-ICE denoising Other ICE tasks Discussion

Outline

1

TV denoising with Iterated Conditional Expectations

2

Other (imaging?) tasks with ICE Deblurring and inverse problems regularized with TV TV-ICE denoising for Poisson noise ICE of a convex functional ICE of a convex set

TV-ICE denoising Other ICE tasks Discussion

The idea of TV-ICE restoration

Definition of the ICE sequence Start with an arbitrary u0 and for all n ∈ N set un+1(x) = 1 Z

• R

un(x)e− Aun−v2+λTV (un)

2σ2

dun(x). computable? convergence un → ˆ uICE? The iterations are easy to deduce from the denoising case! Case where Aδx2 does not depend on x : we have

• wn+1 = un − γA∗(Aun − v)

un+1 = Fwn+1(un) with parameters (γσ2, γλ), where γ = Aδx−2.

TV-ICE denoising Other ICE tasks Discussion

Convergence condition

Assumptions γ = Aδx−2 does not depend on x A1Ω = 1Ω′ Theorem If γ < 2, then (un)n∈N linearly converges to a limit ˆ uICE independent of u0 such that ∀x ∈ Ω, ˆ uICE(x) = Eu|v[u(x) | u(xc) = ˆ uICE(xc)]. But for each γ ≥ 2 there are always cases of non-convergence. deconvolution: if A = k ∗ · with k = 1, then γ = 1/k2. Gaussian blur: γ < 2 ⇔ σA 0.5 pixel → k should be very concentrated. zooming: if A = block-averaging, blocks should have size < 2. → very limited applications!

TV-ICE denoising Other ICE tasks Discussion

4 possible strategies to ensure convergence

• wn+1 = un − γA∗(Aun − v) = (I − γA∗A)un + γA∗v

(1) un+1(x) = Fwn+1(un) (2) 1st strategy: averaging on u. Replace (2) step with un+1 = (1 − r)un + r Fwn+1(un) Observation: r ≤ min(1,

2 γρ(A∗A))

= ⇒ linear convergence. 2nd strategy: averaging on w. Replace (1) step with wn+1 = (1−s)wn +s (un −γA∗(Aun −v)) Observation: s ≤ min(1,

2 γρ(A∗A))

= ⇒ linear convergence. 3rd strategy: set γ free. Replace (1) step with wn+1 = un − τA∗(Aun − v), τ > 0. Observation: τ < 2 = ⇒ linear convergence. 4th strategy: “implicitize”. Replace (1) step with with wn+1 = (I + γA∗A)−1(un + γA∗v). Observation: linear convergence.

TV-ICE denoising Other ICE tasks Discussion

Application to image deblurring

Framework A = k ∗ · ⇒ γ = 1/k2 If k2 > 1/2, the natural strategy applies

• wn+1 = un − γˇ

k ∗ (k ∗ un − v) un+1 = Fwn+1(un) else, the averaging and free-gamma strategies always apply

• wn+1 = (1 − s)wn + s (un − τˇ

k ∗ (k ∗ un − v)) un+1 = (1 − r)un + r Fwn+1(un) Implicit strategy available only for periodic boundary conds:

• wn+1 = (I + γA∗A)−1(un + γA∗v) = F−1

Fu+γ(Fk)∗·F(v) 1+γ|Fk|2

• un+1 = Fwn+1(un).

TV-ICE denoising Other ICE tasks Discussion

Periodic constant blur on a 3.5-ray disc + Gauss. noise with sd. 2

TV-ICE denoising Other ICE tasks Discussion

ROF deblurring

TV-ICE denoising Other ICE tasks Discussion

TV-ICE deblurring by averaging on u (γ ≈ 209)

TV-ICE denoising Other ICE tasks Discussion

TV-ICE deblurring by averaging on w (γ ≈ 209)

TV-ICE denoising Other ICE tasks Discussion

TV-ICE deblurring by free-gamma (γ ≈ 209)

TV-ICE denoising Other ICE tasks Discussion

TV-ICE deblurring by implicit scheme (γ ≈ 209)

TV-ICE denoising Other ICE tasks Discussion 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1 100 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 ROF (primal-dual algorithm) TV-ICE (u-averaging strategy) TV-ICE (w-averaging strategy) TV-ICE (free-gamma strategy) TV-ICE (implicit strategy)

Convergence curves

TV-ICE denoising Other ICE tasks Discussion

Application to zooming (from block-averaging)

Framework: Let z be a zoom factor (z ∈ N∗). A = averaging on z × z-blocks + subsampling by factor z: Au(i, j) = 1 z2

z−1

• k=0

z−1

• l=0

u(zi + k, zj + l). A∗A = averaging on z × z-blocks with no subsampling. Remark: As (I + γA∗A)−1(u + γA∗v) = (I −

γ 1+γ A∗A)u + γ 1+γ A∗v,

3rd and 4th strategies are equivalent when τ =

γ 1+γ < 2.

Algorithm:

• ∀x, wn+1(x) = un(x) +

γ 1+γ (v(x/z) − ¯

un

x)

∀x, un+1(x) = Fwn+1(un).

TV-ICE denoising Other ICE tasks Discussion

zero-order interpolation (4 × 4 zooming)

TV-ICE denoising Other ICE tasks Discussion

bicubic interpolation (4 × 4 zooming)

TV-ICE denoising Other ICE tasks Discussion

ROF (4 × 4 zooming)

TV-ICE denoising Other ICE tasks Discussion

TV-ICE (4 × 4 zooming)

TV-ICE denoising Other ICE tasks Discussion

1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1 100 200 400 600 800 1000 1200 1400 1600 1800 2000 LSE dual ROF primal-dual ROF u-averaged ICE w-averaged ICE free-gamma ICE implicit ICE

λ = 1 and σ = 5 for each algorithm (2 × 2 zooming).

TV-ICE denoising Other ICE tasks Discussion

Outline

1

TV denoising with Iterated Conditional Expectations

2

Other (imaging?) tasks with ICE Deblurring and inverse problems regularized with TV TV-ICE denoising for Poisson noise ICE of a convex functional ICE of a convex set

TV-ICE denoising Other ICE tasks Discussion

TV Poisson denoising

[R. Abergel, C.L., L. Moisan, T. Zeng, SSVM 2015] Poisson noise modelling P(v|u) =

• x∈Ω

u(x)v(x) v(x)! ∝ exp(−u − v log u, 1Ω) So TV denoising posterior p.d.f. is written as P(u|v) = e−u−v log u,1Ω−λTV (u).            ˆ uMAP = arg minuu − v log u, 1Ω + λTV (u) ˆ uLSE(x) =

• (R+)Ω u(x)e−u−v log u,1Ω−λTV (u) du
• (R+)Ω u(x)e−u−v log u,1Ω−λTV (u) du

ˆ uICE = lim un where un+1(x) =

• R+ sv(x)+1e−(s+λ

y∈Nx |s−un(y)|) ds

• R+ sv(x)e−(s+λ

y∈Nx |s−un(y)|) ds

TV-ICE denoising Other ICE tasks Discussion

One ICE iteration: closed formula

Closed form un+1(x) =

• 1≤k≤5 ckI µk,v(x)+1

ak−1,ak

• 1≤k≤5 ckI µk,v(x)

ak−1,ak

, where 0 = a0 ≤ a1 ≤ a2 ≤ a3 ≤ a3 ≤ a5 = +∞ with {a1, a2, a3, a4} = {un(Nx)}, µk = 1 − (6 − 2k)λ, ck = eλ(

k−1

j=1 aj−4 j=k aj) ,

and I µ,p

x,y =

y

x

spe−µsds . It “suffices” to have good computation schemes of upper generalized gamma function and incomplete Gamma function γµ(p, x) = x spe−µsds, Γµ(p, x) = +∞

x

spe−µsds.

TV-ICE denoising Other ICE tasks Discussion

reference

TV-ICE denoising Other ICE tasks Discussion

noisy

TV-ICE denoising Other ICE tasks Discussion

MAP

TV-ICE denoising Other ICE tasks Discussion

ICE

TV-ICE denoising Other ICE tasks Discussion

Convergence + no-staircasing

Theorem The sequence (un) started at u0 = 0 converges linearly to ˆ uICE. No-staircasing result Let x and y ∈ Ω. If ˆ uICE is constant on Nx ∪ Ny ∪ {x, y}, then v(x) = v(y).

TV-ICE denoising Other ICE tasks Discussion

Outline

1

TV denoising with Iterated Conditional Expectations

2

Other (imaging?) tasks with ICE Deblurring and inverse problems regularized with TV TV-ICE denoising for Poisson noise ICE of a convex functional ICE of a convex set

TV-ICE denoising Other ICE tasks Discussion

Definition of ICE for a convex functional J

Let J : Rd → R ∪ {+∞} be a convex, l.s.c. and coercive function with nonempty domain (dom(J) := {J < +∞}). Definition An ICE point of J is a point x ∈ dom(J) such that ∀1 ≤ i ≤ d, xi =

• R xie−J(x)dxi
• R e−J(x)dxi

. The ICE algorithm is a sequence (xn)n∈N started at x0 ∈ dom(J) and such that ∀1 ≤ i ≤ d, ∀n ∈ N, xn+1

i

=

• R xn

i e−J(xn)dxn i

• R e−J(xn)dxn

i

Uniqueness? Existence? Convergence?

TV-ICE denoising Other ICE tasks Discussion

Sufficient conditions for convergence

Theorem

1 If J is C2 with dom(J) = Rd, and if its Hessian H is uniformly

strictly diagonally dominant, then the ICE point exists, is unique and the ICE algorithm converges linearly.

2 If J depends on an image v (so Jv ← J), and if

u → Jv(u) and v → Jv(u) are subdifferentiable; for every 1 ≤ i, j ≤ d, ˇ ui ∈ Rd−1, v ∈ Rd, we have ui → ∂Jv/∂uj and ui → ∂Jv/∂vj nonincreasing; for all u ∈ Rd and c ∈ R, Jv+c(u + c) = Jv(u); for all u and v ∈ Rd, J−v(−u) = Jv(u); for all u ≥ 0 and v ≥ 0 in Rd, Jv(−u) ≥ Jv(u);

then the ICE point exists, is unique and the ICE algorithm converges linearly. Example 1: J(x) = 1

2xTHx − bx ⇒ ICE point=min J and ICE algo

= gradient descent with Jacobi preconditioning ⇒ converges when H is strictly diag. dominant.

TV-ICE denoising Other ICE tasks Discussion

Sufficient conditions for convergence

Theorem

1 If J is C2 with dom(J) = Rd, and if its Hessian H is uniformly

strictly diagonally dominant, then the ICE point exists, is unique and the ICE algorithm converges linearly.

2 If J depends on an image v (so Jv ← J), and if

u → Jv(u) and v → Jv(u) are subdifferentiable; for every 1 ≤ i, j ≤ d, ˇ ui ∈ Rd−1, v ∈ Rd, we have ui → ∂Jv/∂uj and ui → ∂Jv/∂vj nonincreasing; for all u ∈ Rd and c ∈ R, Jv+c(u + c) = Jv(u); for all u and v ∈ Rd, J−v(−u) = Jv(u); for all u ≥ 0 and v ≥ 0 in Rd, Jv(−u) ≥ Jv(u);

then the ICE point exists, is unique and the ICE algorithm converges linearly. Example 2: Jv(u) = Poisson noise(v, u) + TV (u) ⇒ linear convergence.

TV-ICE denoising Other ICE tasks Discussion

Outline

1

TV denoising with Iterated Conditional Expectations

2

Other (imaging?) tasks with ICE Deblurring and inverse problems regularized with TV TV-ICE denoising for Poisson noise ICE of a convex functional ICE of a convex set

TV-ICE denoising Other ICE tasks Discussion

Definition of ICE for a convex set C

Consider the canonical basis (ei)1≤i≤d of Rd (or another basis). Let C ⊂ Rd be a nonempty compact convex set. Definition An ICE point of C (relatively to the basis (ei)) is a point x ∈ C that is the midpoint of [ai(x), bi(x)] for each 1 ≤ i ≤ d, where [ai(x), bi(x)] := C ∩ {x + tei, t ∈ R}. The ICE algorithm is a sequence (xn)n∈N started at x0 ∈ C and such that xn+1 is the middle of [ai(xn), bi(xn)] (1 ≤ i ≤ d). existence of an ICE point by Schauder fixed-point theorem; a rectangle may have infinitely many ICE points; counterpart for center of gravity; translation- but not rotation-invariant.

TV-ICE denoising Other ICE tasks Discussion

Definition An extremal point of C is a point x ∈ C such that ∀1 ≤ i ≤ d, ∀y, z ∈ C \ {x}, y − x, ei · z − x, ei > 0. We have x extremal point ⇒ x ICE point. ⇒ even (strictly) convex sets may have several ICE points (e.g. B((1, 0), 1) ∩ B((0, 1), 1) ⊂ R2). Observation

1 If C is strictly convex and has no extremal point (e.g. if C has

a C2 boundary), then the ICE point is unique.

2 If C has a C2 boundary then the algorithm converges to an

ICE point and the convergence is linear.

TV-ICE denoising Other ICE tasks Discussion

Case of the triangle

Theorem

1 A triangle has (apart from

its extremal vertices) a unique ICE point.

2 The ICE algorithm (not

initialized on an extremal vertex) converges linearly.

3 ICE point found by a simple

geometric construction. Construction of ICE Location of ICE ∀ rotation right isosceles equilateral

TV-ICE denoising Other ICE tasks Discussion

Concluding remarks

framework for some image restoration problems without energy minimization fast approximation of LSE purely primal algorithm: no big theory + any initialization + nice-to-see convergence raises interesting questions even in low dimension.