2.2 Dual Formulation
2.2.1 Discrete algorithm
By considering the minimization problem in a discrete framework we have the following analog to (2.32)
minp
nkdivhp − λ−1d0k2X: |pi,j| ≤ 1, i, j = 1, . . . No
. (2.37)
Deriving with the same arguments as (2.33)-(2.36) and considering pn+1− pn
k = ∇h
divhpn− λ−1d0
− ∇h
divhpn− λ−1d0
pn+1 (2.38) yields a semi-implicit scheme known as Chambolle’s iteration that is given in [9]
p0= 0, pn+1=
pn+ k∇h
divhp − λ−1d0 1 + k
∇h
divhp − λ−1d0
. (2.39)
The iteration converges for k ≤ 1/8 which is proven in [9]. In practice, convergence is observed by setting k equal to 1/4. This gives a rapid convergence up to 50-100 iterations, and then the speed of convergence decays.
Chapter 3
Staircase reducing models
Section 2.1.1 pointed out some of the areas where the ROF model may suffer from an inefficient restoration. This chapter will first give a brief overview of some of the proposals from the literature that may improve these problems.
Then the main focus will be on the staircasing effect. Two models that deal with this effect will be presented: the Fourth-Order model [10, 12, 26] and the TV-Stokes model [38, 34, 25, 21]. The latter reference is the novelty in research literature and it is the main contribution in this thesis. Both models will be given in dual formulation and the next section will give numerical experiments for these models. It should also be noted that these problems are still challenging problems in image processing and remain to be solved.
The ROF model treats flat regions and edges well, but favours piecewise constant functions in the solution, known as staircase effect. This is an adverse effect in images that contain regions with gradual image variations. To pre-vent this over-sharpening there have been some proposals to minimize different regularizations terms.
A more general version of the ROF model may be given as the following minimization problem
infd
Z
Ω
φ(|∇d|) +λ
2(d − d0) dx, (3.1)
where φ : R+ → R+ is a smooth function. The above minimization has the following regularized Euler equation
div
φ0(p|∇d|2+ β) ∇d
|∇d|
− λ(d − d0) = 0, (3.2) where β > 0 is a small parameter to avoid |∇d| = 0. Pure total variation corresponds to φ(g) = g and φ(g) = g2 corresponds to the H1 norm. The question now is how to combine the two norms to benefit each of them.
In [4] they propose
φ0(|∇d|) = |∇d|p−1 (3.3)
and then choose p close to 1.
The smooth function φ in the minimization problem (3.1) will be more in-teresting when endowed with the following conditions
(C1) : φ(0) = φ0(0) = 0 (3.4)
(C2) : lim
s→+∞
φ(s)
s = 0. (3.5)
The first condition (C1) implies that the first term in (3.1) is quadratic near the origin, thus smoothes the flat areas (|∇d| ≈ 0) with the H1 norm. The second condition (C2) implies that the function is sublinear in growth at infinity. This property will enhance the edges since the cost of edges are low. Using the above conditions for the smooth function φ results in an ill-posed problem since the function is obviously non-convex. A simple example would be choosing φ(s) = s2/(1 + s2). However, the non-convexity is irrelevant for the finite-dimensional case as stated in [2] on page 93, and the numerical experiments often give visually good results. This strategy is throughly presented in [23]
where the authors also gave the dual formulation with a weighted total variation term of the non-convex problem.
In [5, 4] they proposed φ(|∇d|) = |∇d|p(|∇d|)to adapt the behaviour of |∇d|, such that p = 1 near the edges and p = 2 in flat regions. In between these regions they use a fractional norm Hp1, p ∈ (1, 2). A simple example would be p(|∇d|) = 1+2|∇d|2 . However, in [5, 4], examples are given in just the one dimensional case.
3.1 Fourth-Order denoising in dual formulation
One of the early approaches to remove staircasing in image denoising, was to introduce higher-order derivatives or norms into the energy. The first approach to combine the TV norm and H1 norm (square norm of the gradient) is to consider the inf-convolution as the authors Chambolle and Lions proposed in [10]. The resulting minimization is equivalent to
inf where is a given threshold to separate the smooth and discontinuous parts of the image.
The same method can be used to minimize the total variation of the gradi-ent, i.e a second-order functional. Some of the earliest variants of higher-order methods for image denoising can be found in the same paper [10], where they decompose the image into two parts, d = d1+ d2, where d1consists of discontin-uous parts and d2consists of the smooth areas of the image, and then minimize the following problem such that each component is assigned to the appropriate norm, i.e d1∈ BV (Ω) and d2∈ H11(Ω) with ∇d2∈ BV (Ω). The results given in [10] are very good for the one-dimensional case, but the model suffers from over-smoothening of the edges in two dimensions.
This is evident, since the problem is how to combine the two different norms, which are dependent on the gradient, and this operation is unstable.
One of the first approaches to minimize only the higher-order derivaives in the energy was proposed by the the authors Lundervold, Lysaker and Tai (LLT) in [26], where they gave the following minimization problem
Z
Ω
(|dxx| + |dyy|) dx dy subject to kd − d0k2≤ σ (3.8) and another functional that is rotational invariant
Z
Ω
q
|dxx|2+ |dxy|2+ |dyy|2+ |dyx|2dx dy subject to kd − d0k2≤ σ. (3.9) The Euler equations for the two above functionals will result into a set of fourth-ordered nonlinear PDEs.
We will further focus on the Fourth-Order model, and derive a dual formu-lation of the following minimization problem
inf regulariza-tion to avoid the non-differentiability. The model shares the same numerical difficulities as the ROF model, and a dual formulation will now be given to overcome the β-regularization and the slow iteration. Similar works can be found in [15].
We can define the regularization term as an analog to the total variation defined in (2.7), hence Thus, we use the right-handed definition of the above equation and deduce
inf where p is the dual variable. By the convexity of d and the concave dual function p, we can use the same argument as in (2.25)-(2.29) to switch the inf and sup in (3.12) Solving for the inner minimization with respect to d yields
d = d0− λ∆p. (3.14)
Thus, by the simple equation above we find our restored image by looking for the projection λ∆p. We deduce a minimization problem for the dual variable by substituting equation (3.14) into (3.13) and setting sup{·} = − inf{−·}, and end up with the following dual problem
infp k∆p − λ−1d0k22: |p| ≤ 1 . (3.15)
The dual minimization problem is solved by considering the optimality con-dition for equation (3.15)
∇2 ∆p − λ−1d0 + αp = 0, (3.16) where either α = 0 and |p| ≤ 1 or α > 0 and |p| = 1, both cases yield
α = |∇2 ∆p − λ−1d0 |. (3.17) Thus, inserting equation (3.17) into (3.16) gives
∇2 ∆p − λ−1d0 + |∆p − λ−1d0|p = 0, (3.18) which may be solved in the same manner as the second-order Chambolle itera-tion
pn+1= pn+ τ ∇2 ∆pn− λ−1d0
1 + τ |∇2(∆pn− λ−1d0)| . (3.19)
3.1.1 Discrete algorithm
The continuous iteration from equation (3.19) is discretized by using the matrix operators from section 1.1.3. Thus, the second-order operator ∇2may have the following discrete definition
∇2hd =
−dBxTBx
−ByTByd
, (3.20)
and the discrete laplacian operator in the same manner
∆hd = −dBTxBx− ByTByd. (3.21) The full discrete version of (3.19) is then given by using the above discrete operators
pn+1=pn+ τ ∇2h ∆hpn− λ−1d0
1 + τ |∇2h(∆hpn− λ−1d0)| . (3.22) The following result gives a condition on τ for fast convergence.
Theorem 2. Let τ ≤ 1/64. Then λ∆hpn converges to solution λ∆hp of the minimization problem given in (3.15) as n → ∞.
A proof of the theorem can be found in [15] where they propose a dual algorithm for the minimization problem given in (3.7).
As an alternative argument for convergence and a better condition for τ , we can consider the following
pn+1= pn+ τ ∇2h ∆hpn− λ−1d0 = (1 + τ ∇2h∆h)pn− τ ∇2h∆hλ−1. (3.23) By denoting the above equation as pn+1 = φ(pn) and requiring φ0(pn) ≤ 1 (1-Lipschitz) we deduce
kτ ∇2h∆hk − k1k ≤ k1 + τ ∇2h∆hk ≤ 1. (3.24) Hence, we get the following condition for τ
τ k∇2h∆hk ≤ 2. (3.25)
To estimate the norm of the operator ∇2h∆h, we note that kBk∞= 2 and the SVD of the operator ∇2h∆h have the following form in x (resp. y) direction
∇2h∆h= −CT
Thus, the following bound in x (resp. y) direction can be estabilished
λmax ∇2h∆h = λmax Σ4x +λmax Σ4xy ≤ kBk42+ kBk42≤ kBk4∞+ kBk4∞= 32, (3.27) giving the optimal condition for τ to be equal to 1/16. Note that the op-erator ∇2h∆h needs to be calculated in each iteration, and this operation is ill-conditioned, i.e the condition number κ(∇2h∆h) is close to 1018.