& %
Department of Mathematics
ROF model on the graph
Japhet Niyobuhungiro and Eric Setterqvist
ROF model on the graph
Japhet Niyobuhungiro∗
Division of Mathematics and Applied Mathematics, Department of Mathematics, Linköping University, SE–581 83 Linköping, Sweden
Department of Mathematics, School of pure and applied sciences, College of Science and Technology, University of Rwanda, P.O. Box 117 Butare, Rwanda
Eric Setterqvist
Division of Mathematics and Applied Mathematics, Department of Mathematics, Linköping University, SE–581 83 Linköping, Sweden
Abstract
In this paper we consider an analogue of the well–known in image processing, Rudin–Osher–Fatemi (ROF) denoising model on a general finite directed and connected graph. We consider the space BV on the graph and show that the unit ball of its dual space can be described as the image of the unit ball of the space
`∞on the graph by a divergence operator. Based on this result, we propose a new fast algorithm to find the exact minimizer for the ROF model. Finally we prove convergence of the algorithm and illustrate its performance on some image denoising test examples.
Keywords: ROF model, Directed graph, L–functional, Image processing, Dual BV, Regularization. 2010 MSC: 46B70, 46E35, 68U10
1. Classical ROF model and ROF model on the graph
Let us suppose that we observed noisy image fob∈L2defined on a square domainΩ= [0, 1]2inR2
fob= f∗+η,
where f∗∈ BV is the original image and η∈L2is noise. Denoising is one of the problems which appear in
image processing: "How to recover the image f∗from the noisy image fob?". Variational methods using the
total variation minimization are often employed to solve this problem. The total variation regularization technique was introduced by Rudin, Osher and Fatemi in [1] and is called the ROF model. It suggests to take as an approximation to the original image f∗ the function fopt,t ∈ BV, which is the exact minimizer
for the L2,1– functional for the couple L2, BV:
L2,1 t, fob; L2, BV = inf g∈BV 1 2kfob−gk 2 L2+tkgkBV , for some t>0, (1) i.e., fopt,t∈BV is such that
L2,1 t, fob; L2, BV = 1 2 fob− fopt,t 2 L2+t fopt,t BV. (2) By BV we denote the space of functions of bounded variation defined by the seminorm
kfkBV(Ω)= Z 1 0 varxf(x, y)dy+ Z 1 0 varyf(x, y)dx, ∗Corresponding author
Email addresses: japhet.niyobuhungiro@liu.se, jniyobuhungiro@ur.ac.rw (Japhet Niyobuhungiro), eric.setterqvist@liu.se(Eric Setterqvist)
where varxf(x, y) = sup 0≤x1≤...≤xn≤1 n−1
∑
j=1 f xj+1, y −f xj, y is a function of y equal to variation of f on the horizontal axis for y fixed, andvaryf(x, y) = sup 0≤y1≤...≤yn≤1 n−1
∑
i=1 |f(x, yi+1) − f(x, yi)|is a function of x equal to variation of f on the vertical axis for x fixed.
However the problem of actual calculation of the function fopt,t(see (1) and (2)) is non-trivial. Standard
approach is connected with discretization of the functional (1), i.e. we divideΩ into N×N square cells and instead of the space L2(Ω) consider its finite dimensional subspace SN which consists of functions
that are constant on each cell.
SN= ( f ∈L2(Ω): f = N
∑
i,j=1 fijχ
ij, ) , (3) whereχ
ij(x, y) = 1 if j−1N ≤x< Nj and i−1N ≤y< Ni; 0 otherwise. (4) Throughout, we will consider our discretization grid as a 2D Cartesian coordinate in (computer) screen space, i.e., with the origin at the top left corner and+y axis pointing down, which is the way matrices are represented on the computer, as illustrated in Figure 1.Figure 1: Illustration of the rectangular grid. The grid contains N2cells refered to as vertices; N(N−1)horizontal edges directed
left right and(N−1)N vertical edges directed up down. Total of 2N(N−1)edges.
As the BV seminorm of a function g∈SNis equal to
kgkBV(SN) = 1 N N
∑
i=1 N−1∑
j=1 gi,j+1−gij + N∑
j=1 N−1∑
i=1 gi+1,j−gij ! ,the discrete analogue of the functional (1) can be written as L2,1 t, fob; L2, BV (5) = inf g∈SN 1 2N2 N
∑
i,j=1 fobij−gij 2 ! + t N N∑
i=1 N−1∑
j=1 gi,j+1−gij + N∑
j=1 N−1∑
i=1 gi+1,j−gij !! .Now we can consider the graph G= (V, E), where the set of vertices V corresponds to cells and the set of edges E corresponds to set of pairs of cells which have common faces (see Figure 1), i.e.
V=vij; i, j=1, 2, . . . , N , and E=EH∪EV, where EH= n eijH= vij, vi,j+1 ; i=1, 2, . . . , N and j=1, 2, . . . , N−1 o are the N(N−1)horizontal edges (that we choose to be directed left right) and
EV =
n
eVij = vij, vi+1,j ; i=1, 2, . . . , N−1 and j=1, 2, . . . , N
o
are the N(N−1)vertical edges (that we choose to be directed top down). So we have that dim V= N2and dim E=2N(N−1). Let us denote by SV and SE the set of real-valued functions on V and E respectively.
Next we consider the analogue of the gradient operator on the graph, i.e., grad : SV −→SE which maps
function f ∈SV to function grad f ∈SE defined as (see Figure 2)
Figure 2: The neighbors of the vertex vij.
(grad f)eijH= f vi,j+1− f vij , i=1, 2, . . . , N and j=1, 2, . . . , N−1;
(grad f)eVij= f vi+1,j− f vij , i=1, 2, . . . , N−1 and j=1, 2, . . . , N.
In this framework, the observed image fob∈SNcan be considered as an element of SV, and the functional
(5) can be written as L2,1 t, fob; L2, BV = inf g∈SV 1 2N2kfob−gk 2 `2(SV)+ t Nkgrad gk`1(SE) , (6)
where kfk`2(SV)=
∑
v∈V (f(v))2 !12 ,khk`1(SE)=∑
e∈E |h(e)| =∑
e∈EH |h(e)| +∑
e∈EV |h(e)|. Of course, exact minimizer of (6) coincides with exact minimizer ofL2,1 s, fob;`2(SV), BV(SV ) = inf g∈SV 1 2kfob−gk 2 `2(SV)+skgrad gk`1(SE) , s=Nt. (7) This leads to the following analogue of the ROF model on a general finite directed and connected graph G= (V, E).
1.1. Formulation of the problem on general graph
Let G = (V, E) be a finite directed and connected graph with N vertices V = {v1, v2, . . . , vN}and M
directed edges E={e1, e2, . . . , eM}, where each edge is determined by a pair of two vertices e= vi, vj for
some i, j∈ {1, 2, . . . , N}. Let SV ={f : f : V−→R}denote the set of of real–valued functions defined on
the vertices and SE ={g : g : E−→R}denote the set of of real–valued functions defined on the edges.
Then the analogue of ROF model on the graph G can be formulated as follows:
Problem 1.1. Suppose that we know function fob∈SV. For given t>0, find exact minimizer of the functional
L2,1 t, fob;`2(SV), BV(SV) = inf g∈BV(SV) 1 2kfob−gk 2 `2(S V)+tkgkBV(SV) , (8) where kfk`2(SV)=
∑
v∈V (f(v))2 !12 ;kfkBV(S V)=kgrad fk`1(SE);khk`1(SE)=∑
e∈E |h(e)|, (9) and operator grad : SV −→SEis defined by the formula(grad f) (e) = f(vj) − f(vi)if e= (vi, vj).
In this paper we will suggest a fast reiterative algorithm that constructs exact minimizer of (8). We would like to note that the case of a general graph could be particularly useful when instead of usually considered rectangular domain we have some manifold, for example map of the Earth is an image on the sphere.
2. Description of the ball of dual space to BV(SV)
It was shown in [2] that exact minimizer for the L2,1-functional
L2,1 t, fob;`2, X = inf g∈X 1 2kfob−gk 2 `2+tkgkX
for the couple `2, X, where space`2is defined by the standard Euclidean normk·k2, X is a Banach space onRn and t is a given positive parameter, is equal to the difference between foband the nearest element
to fobof the ball of radius t>0 of the space X∗. Note that proofs in [2] are also true when X is equipped
with seminorm. So to construct exact minimizer for the L2,1-functional for the couple `2(SV), BV(SV)
we first need to describe the ball of radius t>0 of the space BV∗(SV)with norm defined by
khkBV∗(S V)= sup k f kBV(SV)≤1 hh, fiS V, where hh, fiSV =
∑
v∈V h(v)f(v). (10)Remark 1. Note thatkfkBV(S
V)=0 if and only if f =C for some constant C∈R. Therefore
khkBV∗(S V)= sup k f kBV(SV)≤1 hh, fi , if ∑ v∈V h(v) =0; +∞ , otherwise.
In order to formulate the main result of this section, let us consider divergence operator on the graph, i.e. the operator div : SE →SV defined by
(div g)(v) =
∑
i:(vi,v)∈E
g((vi, v)) −
∑
j:(v,vj)∈Eg (v, vj) . (11)
We can now formulate the result
Theorem 2.1. The unit ball of the space BV∗(SV)is equal to the image of the unit ball of the space`∞(SE)under
the operator div, i.e.,
BBV∗(S V)=div B`∞(SE) .
The proof of Theorem 2.1 is based on two elementary lemmas. The first one is
Lemma 1. The operatordiv : SE−→SV is conjugate to the operator grad : SV −→SE, i.e.
hdiv g, fiS V =hg, grad fiSE, wherehf1, f2iSV = ∑ v∈V f1(v)f2(v)andhg1, g2iSE = ∑ e∈E g1(e)g2(e).
Proof. The expressionhg, grad fiS
E is, by definition of the operator grad, equal to
hgrad f , giS E =
∑
e=(vi,vj)∈E f vj−f(vi) g(e) =∑
e=(vi,vj)∈E f vj g vi, vj−∑
e=(vi,vj)∈E f(vi)g vi, vj =∑
v∈Vi:(v∑
i,v)∈E f(v)g((vi, v)) −∑
v∈Vj:(v,v∑
j)∈E f(v)g v, vj =∑
v∈V f(v) ∑
i:(vi,v)∈E g((vi, v)) −∑
j:(v,vj)∈E g v, vj . The expressionhdiv g, fiSV is, by definition of the operator div, equal to
hf , div giS V =
∑
v∈V f(v) (div g) (v) =∑
v∈V f(v) ∑
i:(vi,v)∈E g((vi, v)) −∑
j:(v,vj)∈E g v, vj . We conclude that hdiv g, fiS V =hg, grad fiSE.For the statement of the second lemma, we will need the following definition
Definition 2.1. Let X be a Banach space and let V be a subspace of X. The annihilator of V denotedAnn(V)is the set of bounded linear functionals that vanish on V. That is the set defined by
Ann(V) ={x∗ ∈X∗: hx∗, vi =0, for all v∈V}, (12) where X∗is the dual space of X.
Lemma 2. Let X be a Banach space with dual space X∗and let V be a finite dimensional subspace of X and Ann(V)
be the annihilator of V. Let x0∈X. Then
inf
v∈Vkx0+vkX =x∗∈Bsup
X∗∩Ann(V)
hx∗, x0i, (13)
Proof. The case x0 ∈ V is obvious. From now on we suppose x0 ∈/ V. Let us take an arbitrary x∗ ∈
BX∗∩Ann(V)and v∈V. Since x∗∈Ann(V)(i.e.,hx∗, vi =0, ∀v∈V) then we have
hx∗, x0i = hx∗, x0+vi ≤ kx∗kX∗kx0+vkX, ∀v∈V.
Since x∗∈ BX∗(i.e.,kx∗kX∗ ≤1) it follows that
hx∗, x0i ≤ kx0+vkX, ∀v∈V.
Therefore since x∗∈ BX∗∩Ann(V)and v∈V were arbitrary, we have that inf
v∈Vkx0+vkX ≥x∗∈Bsup
X∗∩Ann(V)
hx∗, x0i. (14)
In order to prove the reverse inequality, let us consider the space W, which is the algebraic sum between the span of x0and the space V:
W={x0} +V={w∈X : w=λx0+v, v∈V and λ∈R}, (15)
and take v0 ∈ V such that infv∈Vkx0+vkX = kx0+v0kX. The existence of such v0 follows from the
assumption that V is a finite dimensional subspace of X. Without loss of generality we can assume that
kx0+v0kX = 1. Since W is a normed vector space by itself, it is possible to consider its dual space.
Moreover V is a linear subspace of W and then by the Hahn-Banach Theorem, since x0+v0 ∈ W and
x0+v0∈/V there exists a bounded linear functional x∗0∈W∗such that
hx∗0, vi =0 for all v∈V, and hx∗0, x0+v0i =1. (16)
It follows that
x∗0 ∈Ann(V) and hx∗0, x0i =1. (17)
Let us now investigate the action of x∗0on W. Let λ∈ R, v∈V and let w =λx0+v be an element of W.
Then we have hx∗0, wi = hx∗0, λx0+vi = hx∗0, λx0+λv0+v−λv0i =λhx∗0, x0+v0i + hx∗0, v−λv0i =λ, (18) becausex∗ 0, x0+v0 =1 andx∗ 0, v−λv0
=0 since x∗0∈Ann(V)and v−λv0∈V. Let us now describe
the unit ballBW of W. Suppose that w=λx0+v∈ BW where λ6=0. We have that
1≥ kwkX =kλx0+vkX =|λ| x0+ v λ X≥ |λ| kx0+v0kX =|λ|. (19)
Therefore w=λx0+v∈ BWimplies that|λ| ≤1. The W∗- norm of the functional x∗0 is by definition
kx∗0kW∗= sup
w∈BW
hx∗0, wi.
From (18) and (19), it follows that
kx∗0kW∗= sup
w∈BW
λ= sup
|λ|≤1
λ=1. (20)
By invoking the Hahn-Banach theorem, we can extend the functional x∗0:
∃xe0∗∈X∗such that ex∗0|W=x∗0and xe ∗ 0 X∗=kx ∗ 0kW∗=1. (21) From (17) and (21) we conclude that
e
It follows that inf v∈Vkx0+vkX =kx0+v0kX =hx ∗ 0, x0+v0i = D e x∗0, x0+v0 E (22) ≤ sup x∗∈B X∗∩Ann(V) hx∗, x0+v0i = sup x∗∈B X∗∩Ann(V) hx∗, x0i.
Putting (14) and (22) together, we obtain inf
v∈Vkx0+vkX =x∗∈Bsup
X∗∩Ann(V)
hx∗, x0i.
Proof of Theorem 2.1. The operator grad has a kernel which is given by
ker(grad) ={f ∈SV: f =C, for some C∈R}, (23)
and the annihilator (or orthogonal complement) of the ker(grad)is given by Ann(ker(grad)) =
( F∈SV :
∑
v∈V F(v) =0 ) . (24)Since from Lemma 1, div = (grad)∗, where (grad)∗is the conjugate operator of grad, then it is a simple result from linear algebra that
Im(div) =Ann(ker(grad)) and Im(grad) =Ann(ker(div)). (25) From Remark 1, note that if ∑v∈Vh(v) 6= 0 then khkBV∗(S
V) = +∞, so to avoid this case we will
re-strict ourselves on functions h ∈ Ann(ker(grad)). Then from relations (24) and (25), we conclude that BV∗(SV) =Im(div). Therefore for all h∈ BV∗(SV), there exists at least one g∈SEsuch that h=div g.
So let h ∈ BV∗(SV)and fix g0∈ SE such that h =div g0. Now let us consider the infimum of kgk`∞(SE)
over all functions g in the fiber of h∈BV∗(SV)under the operator div.
inf
h=div gkgk`∞(SE)={u∈ker(div)}inf kg0+uk`∞(SE)
Put U=ker(div). Then we have Ann(U) =Im(grad). By applying Lemma 2, we have that inf h=div gkgk`∞(SE)= sup k f∗k`1(SE)≤1, f∗=grad f , f ∈SV hf∗, g0iSE = sup kgrad f k`1(SE)≤1 hgrad f , g0iSE.
From Lemma 1, h=div g0, (9) and (10) follows that
inf h=div gkgk`∞(SE)= sup kgrad f k`1(SE)≤1 hf , div g0iSV = sup kgrad f k`1(SE)≤1 hf , hiS V = sup k f kBV( SV)≤1 hf , hiS V =khkBV∗(S V). It means that khkBV∗(S
V)≤1 if and only if h=div ginf kgk`∞(SE)≤1. (26)
Note that the infimum is attained because ker(div)is a finite dimensional subspace of SE. Therefore there
exists an element
gh∈g0+ker(div), such that kghk`∞(SE)≤1 and div gh=h.
We conclude that BBV∗(S V)=div B`∞(SE) .
2.1. Corollary for the discrete ROF model
As a corollary of Theorem 2.1 we will have the following result (for a reminder of space SN, see
expression (3)).
Corollary 1. Let h ∈ SN. Then khkBV∗(S
N) ≤ 1 if and only if h can be decomposed as h = h1+h2, with
h1, h1∈SNsuch that sup 0≤x≤1 Z x 0 h1(s, y)ds ≤1 and Z 1 0 h1(s, y)ds=0, for all y∈ [0, 1]; and sup 0≤y≤1 Z y 0 h2(x, t)dt ≤1 and Z 1 0 h2(x, t)dt=0, for all x∈ [0, 1].
Proof. Let us take h ∈ SN such thatkhkBV∗(SN) ≤1. By interpretation of the space SN on the graph, we
have from Theorem 2.1 that
khkBV∗(S
N)≤1 if and only if h=div g for some g∈SEsuch that kgk`∞(SE)≤N,
where h vij = (div g) vij =hgeHi,j−1−geijHi+hgei−1,jV −geViji, i, j=1, . . . , N, with the assumption geHi,0=geV0,j=0. It follows that
h=h1+h2,
where functions h1, h2∈SV are defined as follows
h1 vij
=hgei,j−1H −geHiji and h2 vij
=hgeVi−1,j−geViji. (27) From (27) we deduce that
g ei1H = −h1(vi1); g ei2H = − (h1(vi1) +h1(vi2)); .. . geH i,N−1 = − (h1(vi1) +h1(vi2) +. . .+h1(vi,N−1)); g eH iN = − (h1(vi1) +h1(vi2) +. . .+h1(vi,N−1) +h1(viN)). ∀i=1, 2, . . . , N. (28)
From the the first N−1 equations of (28), we have that sup 1≤k≤N−1 g eikH = sup 1≤k≤N−1 k
∑
j=1 h1 vij , ∀i=1, 2, . . . , N. (29) From the the last equations of (28), we deduce that0=geiNH= −
N
∑
j=1
h1 vij , ∀i=1, 2, . . . , N. (30)
From the correspondance between SV and space SN, the expression (29) implies that
sup 1≤k≤N−1 g eikH = sup 1≤k≤N−1 N Z k N 0 h1(s, y)ds = sup 0≤x≤1 N Z x 0 h1(s, y)ds , ∀y∈ [0, 1]. (31) Sincekgk`∞(SE)≤N, it follows from (31) that
sup 0≤x≤1 N Z x 0 h1 (s, y)ds ≤N, ∀y∈ [0, 1].
Therefore sup 0≤x≤1 Z x 0 h1 (s, y)ds ≤1, ∀y∈ [0, 1]. From (30), we have that
N
∑
j=1 h1 vij =N Z NN 0 h1 (s, y)ds=0, ∀y∈ [0, 1]. Therefore Z 1 0 h1 (s, y)ds=0, ∀y∈ [0, 1]. We conclude that sup 0≤x≤1 Z x 0 h1(s, y)ds ≤1 and Z 1 0 h1(s, y)ds=0, for all y∈ [0, 1].Similar arguments for the function h2and the function g acting on the edges eVij will yield that
sup 0≤y≤1 Z y 0 h2 (x, t)dt ≤1 and Z 1 0 h2 (x, t)dt=0, for all x∈ [0, 1].
3. Algorithm for construction of exact minimizer
Below we present a fast reiterative algorithm for the actual construction of the element fopt,t. This
algorithm is a specialized version of an algorithm constructed in [3] for problems concerning a weighted divergence operator on a general finite directed graph.
As was mentioned in the previous section, we know from [2] that the exact minimizer fopt,t for the
L2,1-functional L2,1 t, fob;`2(SV), BV(SV) = inf g∈BV(SV) 1 2kfob−gk 2 `2(SV)+tkgkBV(SV) is given by fopt,t= fob−eh, (32)
where eh is the nearest element to fob, in the metric of`2(SV), in the ball tBBV∗(S
V), i.e. inf h∈tBBV∗ (SV) kfob−hk`2(SV)= fob−eh `2(S V) . (33)
Next, from Theorem 2.1 we know that tBBV∗(S
V) = t div
B`∞(SE)
for t > 0. The proposed algorithm constructs eh through a sequence of elements gn ∈tB`∞(SE)such that
3.1. Algorithm
The algorithm consists of several steps outlined in this section.
Let fob ∈ SV on G= (V={v1, . . . , vN}, E={e1, . . . , eM}), a regularization parameter t and a maximum
number of iterations Niterbe given. Set
ek= vi, vj
∈E, k=1, 2, . . . , M; for some i, j∈ {1, 2, . . . , N}. Define the operator T as follows:
T=TMTM−1TM−2. . . T2T1,
where for k=1, 2, . . . , M, Tk: tB`∞(SE)−→tB`∞(SE)is defined as follows:
(Tkg) (e) = Kg(ek) ,if Kg(ek) ∈ [−t,+t]; −t ,if Kg(ek) < −t; +t ,if Kg(ek) > +t. ,if e=ek; g(e) ,if e6=ek. where Kg(ek) = fob(vj) − div\ekg (vj)− fob(vi) − div\ekg (vi) 2 . and div\ekg(vi) = (div g) (vi) +g(ek); div\ekg vj= (div g) vj−g(ek); and div\ekg(v`) = (div g) (v`), ∀` 6=i, j.
With these notions at our disposal, we now outline the steps of the algorithm.
Step 1. Take g0=0, or choose any g0∈tB`∞(SE)
Step 2. Calculate g = Tg0. i.e., calculate (Tg0) (ek) for k = 1, 2, . . . , M. If g = g0 then take eh = div(g0),
otherwise go to Step 3.
Step 3. Put g0=g and go to Step 2.
We continue this process applying the operator T to the new element g∈tB`∞(SE)generating the sequence
of elements g0, g1=Tg0, g2=Tg1, . . . , gn=Tgn−1with gn ∈tB`∞(SE), n=0, 1, 2, . . . until the maximum
number of iterations Niteris reached. In the next Section we will show that
div(gn) −→eh as n→ +∞ in the metric of`2(SV).
To prove the convergence of the proposed algorithm we need the following proposition
Proposition 3.1. Let eh be the minimizer defined by (33). The operator T is continuous and satisfies the following two conditions
(a) For any g∈tB`∞(SE), div g=eh if and only if Tg=g; (b) For any g∈tB`∞(SE), if div g6=eh and
kfob−div(Tg)k`2(SV)<kfob−div gk`2(SV). (34) Proof. It is clear from definition that small changes of g∈tB`∞(SE)leads to small changes of any operator
Tkand therefore the operator T is continuous as a product of continuous operators.
We now prove (a). Assume first that g∈ tB`∞(SE)is such that div g =eh. Let ek = (vi, vj) ∈E. We note that g(ek)appears only in the following two terms ofkfob−div gk2`2(S
V): fob(vj) − (div g)(vj)2+ (fob(vi) − (div g)(vi))2= (35) h fob(vj) − (div\ekg)(vj) −g(ek) i2 +hfob(vi) − (div\ekg)(vi) +g(ek) i2 .
Since div g=eh, g(ek)in particular must minimize (35) in the interval[−t, t]. We note that by the Jensen’s inequality ω(g(ek)) = h fob(vj) − (div\ekg)(vj) −g(ek) i2 +hfob(vi) − (div\ekg)(vi) +g(ek) i2 ≥ 2 h fob(vj) − (div\ekg)(vj) −g(ek) i +hfob(vi) − (div\ekg)(vi) +g(ek) i 2 2 = 2 h fob(vj) − (div\ekg)(vj) i +hfob(vi) − (div\ekg)(vi) i 2 2
with equality if and only if
fob(vj) − (div\ekg)(vj) −g(ek) = fob(vi) − (div\ekg)(vi) +g(ek), i.e. when g(ek) = h fob(vj) − (div\ekg)(vj) i −hfob(vi) − (div\ekg)(vi) i 2 = (Kg)(ek).
Moreover, ω(x)is strictly convex and therefore strictly decreasing for x< (Kg) (ek)and strictly increasing
for x > (Kg)(ek). So the minimal value of ω(x) on the interval [−t, t] is only attained at i) the point
(Kg)(ek) if (Kg)(ek) ∈ [−t, t], ii) the point −t if (Kg)(ek) < −t and iii) the point t if(Kg)(ek) > t. The
assumption div g=eh then implies that g(ek)must be the nearest point in the interval[−t, t]to(Kg)(ek), i.e. Tkg(ek) =g(ek). Since ek ∈E was arbitrary it follows that Tkg(ek) =g(ek), k=1, ..., M which implies
Tkg=g, k=1, ..., M and therefore Tg=g.
We now show the other direction of (a). Let us assume that g ∈ tB`∞(SE) and Tg = g. Then for
any edge e ∈ E, g(e) coincides with the point of the interval [−t, t] which is nearest to (Kg)(e). As
kfob−div gk`2(SV)is a convex function of g on tB`∞(SE), to show this direction it is enough to show that g
minimizeskfob−div gk`2(S
V)locally, i.e. it is enough to show that for some small ε>0 we have
kfob−div gk`2(SV)= inf
f ∈Rε
kfob−div fk`2(SV),
where Rε is the tubular set given by
Rε=
n
f ∈tB`∞(SE): kg−fk`∞(SE)≤ε
o . Note that for any f ∈Rεand e∈E we have
f(e) ∈ [−t, t] ∩ [g(e) −ε, g(e) +ε].
The set Rεis a closed convex subset of SE, therefore there exists a function fε ∈Rε such that kfob−div fεk`2(S
V)= f ∈Rinf ε
kfob−div fk`2(S V)
so we need to prove that
kfob−div gk`2(SV)=kfob−div fεk`2(SV).
We first note that it follows from the necessity direction proved above that for any edge e∈ E, fε(e)will
coincide with the point of the interval[−t, t] ∩ [g(e) −ε, g(e) +ε]which is nearest to(K fε)(e).
Let us now decompose the set of edges E into two parts, the first part which we denote byΩgconsists
of the edges for which(Kg)(e)does not belong to the interval[−t, t], i.e. Ωg={e∈E :(Kg)(e)∈ [−/ t, t]}.
As g(e)is the nearest point in the interval[−t, t]to(Kg)(e)we have g(e) =
−t ,if Kg(e) < −t;
If the number ε > 0 is small enough, it follows fromkg− fk`∞(SE) ≤ εthat on e ∈ Ωg where we have
(Kg)(e) < −t we will also have(K fε)(e) < −t and therefore fε(e) = −t= g(e). Analogously, if ε >0 is
small enough, on e∈ Ωgwhere(Kg)(e) >t we will have(K fε)(e) > t and therefore fε(e) = t=g(e). So
we have proved that
fε(e) =g(e)for all e∈Ωg. (36)
The next step is to consider the graph G0 = (V, E\Ωg), i.e. the graph G with edges in Ωg removed.
G0 is the union of several connected components (Vk, Ek), k = 1, ...,` (we have V1∪...∪V` = V and
E1∪...∪E`=E\Ωg). Note that some of the graphs(Vk, Ek)can consist just of one single vertex. For these
graphs there is nothing to prove as Ek=∅.
Let us now consider a subgraph(Vk, Ek)where Ek 6=∅. On each e∈ Ek we have(Kg)(e) ∈ [−t, t]and
therefore g(e) = (Kg)(e), i.e. if e= (vi, vj)then
g(e) = (Kg)(e) = h fob(vj) − (div\eg)(vj) i −hfob(vi) − (div\eg)(vi) i 2 , or equivalently fob(vj) − (div g)(vj) = fob(vi) − (div g)(vi).
(Here it is important to note that operators K, div and div\e are considered in the orignal setting of
G= (V, E).) Therefore, for all v∈Vkthe values of
fob(v) − (div g)(v)
are equal. So we have
∑
v∈Vk [fob(v) − (div g)(v)]2=|Vk| ∑v∈Vk[fob(v) − (div g)(v)] |Vk| !2 . (37)For fεwe can with Jensen’s inequality derive the corresponding inequality
∑
v∈Vk [fob(v) − (div fε)(v)] 2≥ | Vk| ∑v∈Vk[fob(v) − (div fε)(v)] |Vk| !2 . (38)Now, note that flow on edges in Ekare cancelled in the sums
∑
v∈Vk
[fob(v) − (div g)(v)] and
∑
v∈Vk[fob(v) − (div fε)(v)]. (39)
Therefore, only flows on edges inΩgremain in these sums. It then follows from (36) that
∑
v∈Vk
[fob(v) − (div fε)(v)] =
∑
v∈Vk
[fob(v) − (div g)(v)].
Therefore, taking into account (37) and (38), we obtain
∑
v∈Vk [fob(v) − (div fε)(v)] 2≥∑
v∈Vk [fob(v) − (div g)(v)]2.Summing over all Vkthen gives
kfob−div fεk 2 `2(S V)≥ kfob−div gk 2 `2(S V).
Hence, from the definition of fεwe must have kfob−div fεk
2
`2(SV)=kfob−div gk2`2(SV).
To prove (b) note that for ∀g ∈tB`∞(SE)the operators Tk, k =1, ..., M, by definition (see also necessity
part of (a)) satisfy
kfob−div(Tkg)k`2(SV)≤ kfob−div gk`2(SV)
with equality if and only if Tkg(ek) =g(ek). This implies for the operator T=TMTM−1...T1that
kfob−div(Tg)k`2(S
V)≤ kfob−div gk`2(SV)
with equality if and only if Tg= g which in turn by condition (a) is equivalent to div g =eh. Hence for any g∈tB`∞(SE), if div g6=eh then
kfob−div(Tg)k`2(S
V)<kfob−div gk`2(SV),
which establishes (b) and completes the proof.
3.2. Convergence of the algorithm
With Proposition 3.1 at hand, we can now prove convergence of the Algorithm.
Theorem 3.1. Let eh be the minimizer defined by (33), g∈tB`∞(SE)and let T be the operator constructed in Section 3.1. Then
div(Tng) −→eh as n→ +∞ in the metric of`2(SV). (40)
Proof. Let T be the operator constructed in Section 3.1. From Proposition 3.1, T : tB`∞(SE)−→tB`∞(SE)is continuous and satisfies the conditions(a)and(b). Then
Tng∈tB`∞(SE), ∀g∈tB`∞(SE); n=0, 1, 2, . . . . (41) Since T satisfies conditions (a) and (b), the sequence of numbers (fob−div(Tng))n∈N
`2(S
V) is
mono-tonically decreasing and is bounded below by inf
g∈tB`∞(SE)kfob−div gk`2(SV). Therefore it converges. Let us
now consider the sequence(Tng)n∈N⊂tB`∞(SE). Since tB`∞(SE)is a compact set, then(Tng)n∈Ncontains a convergent subsequence in tB`∞(SE), say(Tnkg)
k∈N:
lim
k→∞(T
nkg) =g
h∈tB`∞(SE). (42)
Since T, div andk·k`2(S
V)are continuous operators, we have
kfob−div(Tgh)k`2(S V)= fob −div T lim k→∞(T nkg) `2(S V) = lim k→∞kfob−div(T(T nkg))k `2(S V) = lim k→∞ fob−div Tnk+1g `2(SV).
From Proposition 3.1 and continuity of div andk·k`2(SV)it follows that
kfob−div(Tgh)k`2(S V)=k→∞lim fob−div Tnk+1g `2(S V) ≥ lim k→∞kfob−div(T nk+1g)k `2(S V) = fob−div lim k→∞(T nk+1g) `2(S V) =kfob−div ghk`2(S V).
Therefore by Proposition 3.1 we conclude that div gh=eh. By continuity of div andk·k`2(S
V)we have that lim k→∞kfob−div(T nkg)k `2(S V)=kfob−div ghk`2(SV).
Now, since the subsequencekfob−div(Tnkg)k`2(S V)
k∈Nmust converge to the same limit as
kfob−div(Tng)k`2(SV) n∈N, we conclude that lim n→∞kfob−div(T ng)k `2(S V)=kfob−div ghk`2(SV)= fob−eh `2(S V) .
Therefore, as eh is the unique nearest element to fobin tBBV∗(S
V)we conclude that
div(Tng) −→div gh=h as ne → +∞ in the norm of`2(SV).
4. Tikhonov functional
In this section, we consider the Tikhonov functional in order to use it for comparaison purposes in the next section. Let us consider the Tikhonov functional for the couple of Sobolev spaces L2(Ω), ˙W1,2(Ω) defined on a square domainΩ= [0, 1]2inR2. It is given by:
L2,2 t, fob; L2(Ω), ˙W1,2(Ω) = inf g∈ ˙W1,2(Ω) 1 2kfob−gk 2 L2(Ω)+tkgk2W˙1,2(Ω) , (43) where kfkW˙1,2(Ω) = ∂ ∂xf(x) L2(Ω) = ∂ ∂x1 f(x1, x2) 2 L2 + ∂ ∂x2 f(x1, x2) 2 L2 !1/2 , (44) and kfkL2(Ω) = Z 1 0 Z 1 0 (f(x1, x2)) 2 dx1dx2 1/2 . (45)
We need to find the element fopt,t∈W˙1,2(Ω), exact minimizer for (43). That is fopt,t∈W˙1,2(Ω)such that
L2,2 t, fob; L2(Ω), ˙W1,2(Ω) = 1 2 fob−fopt,t 2 L2(Ω)+t fopt,t 2 ˙ W1,2(Ω). (46)
The exact minimizer fopt,t ∈ W˙ 1,2 for (43) can be found directly by solving the corresponding optimality
condition. However functions from SN don’t have derivative (44). So we need to look for an approximate
smooth function for which (44) exists. We consider a finite dimensional subspace of ˙W1,2(Ω) which consists of functions that can be represented as trigonometric polynomials of degree≤N. To this end, we will carry the processing in the Fourier transform domain where the image is represented as a weighted sum of the complex exponentials. Afterwards, we apply the inverse Fourier transform to return to the spatial domain. We will consider
fob(n1, n2) = 1 N2 N−1
∑
k=0 N−1∑
`=0 cfob(k,`)e
i2π N (kn1+`n2), n 1, n2=0, 1, . . . , N−1, (47) and g(n1, n2) = 1 N2 N−1∑
k=0 N−1∑
`=0 b g(k,`)e
i2πN(kn1+`n2), n 1, n2=0, 1, . . . , N−1, (48) where c fob(k,`) = N−1∑
n1=0 N−1∑
n2=0 fob(n1, n2)e
− i2π N (kn1+`n2), k,` =0, 1, . . . , N−1, (49) and b g(k,`) = N−1∑
n1=0 N−1∑
n2=0 g(n1, n2)e
− i2π N(kn1+`n2), k,` =0, 1, . . . , N−1. (50)Equations (49)–(50) correspond to the discrete Fourier transforms of the images defined over the discrete finite 2D grid of size N×N and equations (47)–(48) correspond to their inverses. We use N samples in both directions and in both spatial and frequency domains. Note that the image is periodized along both dimensions with period N. In this formulations, the functional (43) becomes
inf b g(k,`) 1 2 N−1
∑
k,`=0 cfob(k,`) −gb(k,`) 2 +t "N−1∑
k,`=0 i 2π Nk 2 |bg(k,`)|2+ N−1∑
k,`=0 i 2π N` 2 |gb(k,`)|2 #! . So we need to minimize each term inside the sum, i.e.,inf b g(k,`) 1 2 cfob(k,`) −gb(k,`) 2 +t 2π Nk 2 + 2π N` 2! |bg(k,`)|2 ! , k,` =0, 1, . . . , N−1. (51) We will represent the complex numbers cfob(k,`)andgb(k,`)in their polar form, i.e.,
c
fob(k,`) =rk,`
e
iϕk,` andbg(k,`) =erk,`e
iϕk,`. (52)
Note that rk,` and ϕk,`, k,` = 0, 1, . . . , N−1 are known because fob is given. We only need to solve for
erk,`’s. Putting (52) in (51), we need to solve inf erk,` 1 2(rk,`−erk,`) 2+ 2π N 2 tk2+ `2er2k,` ! , k,` =0, 1, . . . , N−1. Taking derivative with respect toerk,`, we see that the infimum is attained aterk,`which satisfies
− (rk,`−erk,`) + 2π N 2 tk2+ `2· 2erk,`=0, k,` =0, 1, . . . , N−1, we have erk,`= rk,` 1+2 2πN2 t(k2+ `2), k,` =0, 1, . . . , N−1.
Therefore the exact minimizer fopt,tis given by the formula
fopt,t(n1, n2) = 1 N2 N−1
∑
k=0 N−1∑
`=0 [fopt,t(k,`)e
i2π N (kn1+`n2), n 1, n2=0, 1, . . . , N−1, (53) where [fopt,t(k,`) = rk,` 1+2 2πN2 t(k2+ `2)e
iϕk,`. The matlab routines f f t2 and i f f t2 are employed for a fast computation of (53).5. Experimental results
In this section we will illustrate on some test images the performance of ROF model using our algorithm and we will compare with result obtained by using the Tikhonov functional described in the previous section. All the test observed images used have width 512 and height 512 pixels and the corresponding original image is provided for comparison purposes. In view of our analogue of ROF model on the graph (see (7)), note that the regularization parameter we use for ROF model is s=Nt, where N =512 and t is the regularization parameter appearing in the original model (see (1)). The original image used in Figures 3–4 is an 8–bit gray–scale image with intensity values ranging from 0 (black) to 255 (white). The original image used in Figures 5–6 is an artificial image with 4 parts of intensity values 180, 120, 90 and 50. The original image used in Figures 7–8 is an artificial binary image of intensity values 150 and 100. In all the cases the noisy image is the corrupted version of the original image where a Gaussian additive noise of
standard deviation ε with different values of ε is added to the original image. i.e, if we denote the original image by f∗, then the observed noisy image fobis defined by
fob= f∗+ε×randn(512, 512),
where randn(N, N) is a function returning a N–by–N matrix containing pseudorandom values drawn from the standard normal distribution as buit in Matlab.
In Figures 4, 6 and 8, we have highlighted in white the parts of the reconstructed image fopt,t where
f∗i,j−fopt,ti,j
≥α, i, j=1, 2, . . . , N, for different values of α in both reconstructions.
Figure 3: Lenna, the original is a 512×512 image with intensity values ranging from 0 to 255. Top left: Original image f∗; Top
right: Noisy image fob=f∗+ε×randn(512,512), i.e image with Gaussian additive noise of standard deviation ε=23.6. Bottom left:
Figure 4: Top: Residual image, i.e., the difference between the observed image and the reconstructed image fob−fopt,t. Top
left: Residual from Tikhonov reconstruction. Top right: Residual from ROF reconstruction. Bottom: The white shows, matrix componentwise, parts of the recontructed image fopt,twhere
f∗i,j−fopt,ti,j
≥35, i, j=1, 2, . . . , N. Left: in Tikhonov reconstruction. Right: in ROF reconstruction.
Figure 5: Geometric features with stairs, the original is a 512×512 image with 4 regions of intensity values 180, 120, 90 and 50. Top left: Original image f∗; Top right: Noisy image fob = f∗+ε×randn(512,512), i.e image with Gaussian additive noise of standard
Figure 6: Geometric features with stairs. Top: Residual image, i.e., the difference between the observed image and the reconstructed image fob−fopt,t. Top left: Residual from Tikhonov reconstruction. Top right: Residual from ROF reconstruction. Bottom: The
white shows, matrix componentwise, parts of the recontructed image fopt,twhere
f∗i,j−fopt,ti,j
≥20, i, j=1, 2, . . . , N. Left: in Tikhonov reconstruction. Right: in ROF reconstruction.
Figure 7: Geometric features with a cusp, the original is a 512×512 image with regions of intensity values 150 and 100. Top left: Original image f∗; Top right: Noisy image fob=f∗+ε×randn(512,512), i.e image with Gaussian additive noise of standard deviation ε=30.85. Bottom left: Tikhonov reconstruction using t=4.37. Bottom right: ROF reconstruction using s=36.5.
Figure 8: Geometric features with a cusp. Top: Residual image, i.e., the difference between the observed image and the reconstructed image fob−fopt,t. Top left: Residual from Tikhonov reconstruction. Top right: Residual from ROF reconstruction. Bottom: The
white shows, matrix componentwise, parts of the recontructed image fopt,twhere
f∗i,j−fopt,ti,j
≥20, i, j=1, 2, . . . , N. Left: in Tikhonov reconstruction. Right: in ROF reconstruction.
6. Final remarks and discussions
Since its appearance in 1992, the ROF model has received a large amount of popularity for its effeciency in regularizing images without smoothing the boundaries, and it has since been applied to a multitude of other imaging problems (see for example the book [4]). Some of the earlier works for minimizing the total variation based on dual formulation include [5], where the methods presented is based on removal of some of the singularity caused by the non-differentiability of the quantity grad f appearing in the regularization term. In 2004, Chambolle (see [6]) provided an algorithm related to [5] for minimizing the total variation of an image and proved its convergence. In the last few years image decomposition models into a piecewise– smooth and oscillating components that usually people refer to as cartoon and textures (or textures +
noise) respectively, have received interest in the image processing community. For example fopt,t ∈ BV
satisfying (2) is such that
fob= fob− fopt,t+ fopt,t.
This is the decomposition of fob into the piecewise–smooth component fopt,t ∈ BV and the component
fob−fopt,t ∈ L2 wich contains textures and noise. The original theoretical model for such an image
decomposition was introduced in 2002 by Yves Meyer in [7] by using the total variation to model the piecewise–smooth component and an appropriate dual space named G which is the Banach space com-posed of the distributions f which can be written f =∂1g1+∂2g2=div(g), with g1and g2in L∞(Ω)with
the following norm
kfkG=inf
kgkL∞(Ω;R2): f =div(g), kgkL∞(Ω;R2)=ess supx∈Ω
q
|g1(x)|2+|g2(x)|2
,
to model the oscillating component. Some of the most known works proposed in the literature for numeri-cally solving the Meyer’s model or its variants include for instance the works of [8] who proposed a model that splits the image into three components, a geometrical components modeled by the total variation, a texture component modeled by a negative Sobolev norm and a noise component modeled by a negative Besov norm; the works of [6, 5, 9, 10] that proposed an efficient projection algorithm to minimize the total variation. In 2006 [11] and [12] adapted Chambolle’s algorithm in cases of presence of an operator like a convolution kernel for example. Recently, [13] showed that the Bregman iteration is a very efficient and fast way to solve TV problems among other L1–regularized optimization problems and proposed a split Bregman method which they applied to the Rudin–Osher–Fatemi functional for image denoising and [14] designed an algorithm by using the Split Bregman iterations and the duality used by Chambolle to find the minimizer of a functional based on Meyer G–norm. Other works based on the Meyer’s G–norm include for example [15] and [16].
Acknowledgement. The authors are very grateful to Professor Natan Kruglyak for his invaluable suggestions and
discussions that made a considerable improvement to this work.
References
[1] L. I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica, North-Holland D (60) (1992) 259–268.
[2] N. Kruglyak, J. Niyobuhungiro, Characterization of Optimal Decompositions in Real Interpolation, submitted.
[3] N. Kruglyak, E. Setterqvist, Weighted Divergence on the Graph.
[4] T. F. Chan, J. Shen, Image Processing and Analysis. Variational, PDE, Wavelet, and Stochastic Methods, Siam, 2005.
[5] T. F. Chan, G. H. Golub, P. Mulet, A nonlinear primal–dual method for total variation–based image restoration, Siam J. Sc. Comput. 20 (6) (1999) 1964–1977.
[6] A. Chambolle, An algorithm for total variation minimization and applications, Journal of mathemati-cal imaging and vision 20 (2004) 89–97.
[7] Y. Meyer, Oscillating Patterns in Image Processing and Nonlinear Evolution Equations, University Lecture Series, Vol. 22, AMS Providence, 2002.
[8] J. F. Aujol, A. Chambolle, Dual norms and image decomposition, International journal of computer vision 63 (1) (2005) 85–104.
[9] J. F. Aujol, G. Aubert, L. Blanc-Feraud, A. Chambolle, Image decomposition into a bounded variation component and an oscillating component, Journal of mathematical imaging and vision 22 (2005) 71– 88.
[10] A. Chambolle, J. Darbon, Algorithms for total variation minimization, Mathematical model for multi-channel image processing, Beijing 7 (2006) 7–8.
[11] J. Gilles, Noisy image decomposition: a new structure, texture and noise model based on local adap-tivity.
[12] J. F. Aujol, G. Gilboa, T. Chan, S. Osher, Structure-texture image decomposition-modeling, algorithms, and parameter selection, International Journal of Computer Vision 67 (1).
[13] T. Goldstein, S. Osher, The split bregman method for L1–regularized problems, Siam J. Imaging Sci-ences 2 (2) (2009) 323–343.
[14] J. Gilles, S. Osher, Bregman implementation of Meyer’s G–norm for cartoon + textures decomposition. [15] D. M. Strong, J. F. Aujol, T. F. Chan, Scale recognition, regularization parameter selection and Meyer’s
G– norm in total variation regularization, Multiscale Model. Simul 5 (1) (2006) 273–303.
[16] S. Osher, O. Scherzer, G–norm properties of bounded variation regularizaton, Comm. Math. Sci 2 (2) (2004) 237–254.
