第1页 / 共28页
第2页 / 共28页
第3页 / 共28页
第4页 / 共28页
第5页 / 共28页
第6页 / 共28页
第7页 / 共28页
第8页 / 共28页
A nonlinear inverse scale space method for a convex multiplicati....pdf
SIAM J. IMAGING SCIENCES
Vol. 1, No. 3, pp. 294–321
c 2008 Society for Industrial and Applied Mathematics
A Nonlinear Inverse Scale Space Method
∗
for a Convex Multiplicative Noise Model
†
Jianing Shi
‡
and Stanley Osher
Abstract. We are motivated by a recently developed nonlinear inverse scale space method for image denoising
[M. Burger, G. Gilboa, S. Osher, and J. Xu, Commun. Math. Sci., 4 (2006), pp. 179–212; M. Burger,
S. Osher, J. Xu, and G. Gilboa, in Variational, Geometric, and Level Set Methods in Computer Vi-
sion, Lecture Notes in Comput. Sci. 3752, Springer, Berlin, 2005, pp. 25–36], whereby noise can be
removed with minimal degradation. The additive noise model has been studied extensively, using the
Rudin–Osher–Fatemi model [L. I. Rudin, S. Osher, and E. Fatemi, Phys. D, 60 (1992), pp. 259–268],
an iterative regularization method [S. Osher, M. Burger, D. Goldfarb, J. Xu, and W. Yin, Multi-
scale Model. Simul., 4 (2005), pp. 460–489], and the inverse scale space flow [M. Burger, G. Gilboa,
S. Osher, and J. Xu, Commun. Math. Sci., 4 (2006), pp. 179–212; M. Burger, S. Osher, J. Xu, and
G. Gilboa, in Variational, Geometric, and Level Set Methods in Computer Vision, Lecture Notes
in Comput. Sci. 3752, Springer, Berlin, 2005, pp. 25–36]. However, the multiplicative noise model
has not yet been studied thoroughly. Earlier total variation models for the multiplicative noise
cannot easily be extended to the inverse scale space, due to the lack of global convexity. In this
paper, we review existing multiplicative models and present a new total variation framework for the
multiplicative noise model, which is globally strictly convex. We extend this convex model to the
nonlinear inverse scale space flow and its corresponding relaxed inverse scale space flow. We demon-
strate the convergence of the flow for the multiplicative noise model, as well as its regularization
effect and its relation to the Bregman distance. We investigate the properties of the flow and study
the dependence on flow parameters. The numerical results show an excellent denoising effect and
significant improvement over earlier multiplicative models.
Key words.
inverse scale space, total variation, multiplicative noise, denoising, Bregman distance
AMS subject classifications. 35-xx, 65-xx
DOI. 10.1137/070689954
1. Introduction. Image denoising is an important problem of interest to the mathemat-
ical community that has wide application in fields ranging from computer vision to medical
imaging. A variety of methods have been proposed over the last decades, including traditional
filtering, wavelets, stochastic approaches, and PDE-based variational methods. We refer the
reader to [9] for a review of various methods. The additive noise model has been extensively
studied, using the Rudin–Osher–Fatemi (ROF) model [23], an iterative regularization method
[21], and the inverse scale space (ISS) flow [5, 6]. However, the multiplicative noise has not yet
been studied thoroughly. In this paper, we obtain a new convex multiplicative noise model
and extend it to the nonlinear ISS.
∗
Received by the editors April 30, 2007; accepted for publication (in revised form) May 30, 2008; published
electronically September 4, 2008. This research was supported by NSF grants DMS-0312222 and ACI-0321917 and
by NIH G54 RR021813.
http://www.siam.org/journals/siims/1-3/68995.html
†
Department of Biomedical Engineering, Columbia University, 351 Engineering Terrace, MC 8904, 530 W. 120th
St., New York, NY 10027 (js2615@columbia.edu).
‡
Department of Mathematics, UCLA, 520 Portola Plaza, Los Angeles, CA 90095 (sjo@math.ucla.edu).
294
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 12/01/12 to 125.71.231.223. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
INVERSE SCALE SPACE FOR MULTIPLICATIVE NOISE
295
In the additive noise model, one sought to recover the signal u which was corrupted by
additive noise v; i.e., f = u + v. The ROF model [23] introduced total variational minimization
to image processing and defined the solution as follows:
(1.1)
u = arg min
u∈BV (Ω)
|u|BV +
f − u2
L2
λ
2
Ω
variation on Ω, equipped with the BV seminorm which is formally given by |u|BV =
for some scale parameter λ > 0, where BV (Ω) denotes the space of functions with bounded
|∇u|,
also referred to as the total variation (TV) of u. Successful implementations of this mini-
mization problem include the Euler–Lagrange equation by gradient descent on the original
minimization problem [23], second-order cone programming [15], duality [8], and an extremely
fast method based on graph cuts [11]. The choice of the positive scale constant λ in the ROF
model is important in that large λ corresponds to a small amount of noise removal, while
small λ can result in oversmoothing the image. The ROF model was extended to an iterative
regularization method based on the Bregman distance in [21], motivated by Meyer’s analysis
in [19], where he defined texture, “highly oscillatory patterns in an image,” as elements of
the dual space of BV (Ω). In order to preserve the texture information Meyer suggested a
modified variational problem using the space (G,·∗), the dual space for the closure of the
Schwartz space in BV (Ω), equipped with the ∗ (see [19]), as
{|u|BV + λf − u∗}.
(1.2)
u = arg min
u∈BV (Ω)
Meyer also used his analysis to quantify and explain the observed loss of contrast in the ROF
model.
The iterative regularization model, introduced in [21], improved the quality of regularized
solutions in terms of texture preservation as well as the signal-to-noise ratio (SNR). The
nonlinear ISS method devised in [5, 6] formulates the iterative regularization method as a
continuous time flow and provides a better temporal resolution and a better stopping criterion
and is usually much faster to implement than the iterative regularization. Both methods have
solutions which start from an initial condition u0(x) = 0 (assuming
f = 0), come close to the
true u, and then approach the noisy image f . The idea behind these two methods is that larger
scales converge faster than smaller ones, where scale can be precisely defined using Meyer’s
∗. The ISS method has proved to yield better denoising results than standard ROF. In fact
it yields improved denoising results among PDE-based methods [4, 5, 6, 9, 21, 23].
Multiplicative noise has not been as well studied. We consider the problem of seeking the
true signal u from a noisy signal f , corrupted by multiplicative noise η. This arises in medical
imaging, e.g., magnetic field inhomogeneity in MRI [14, 20], speckle noise in ultrasound [28],
and speckle noise in synthetic aperture radar (SAR) images [2, 27]. Denoising of speckled
ultrasound images has been tackled, e.g., in [1], and more general noise models, which are
correlated with the signal amplitude (such as multiplicative noise), have been studied in
[24]. The first total variation approach to solving the multiplicative model was presented in
[22], which used a constrained optimization approach with two Lagrange multipliers. The
authors also considered blurry and noisy images. However, their fitting term is nonconvex,
in general, which leads to difficulties in using the iterative regularization or the ISS method.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 12/01/12 to 125.71.231.223. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
296
JIANING SHI AND STANLEY OSHER
In [2] a multiplicative model was introduced involving BV , with a fitting term derived from
a maximum a posteriori (MAP) estimation. However, the fitting team was convex only for
u∈ (0,2f ). Some interesting analysis for the convergence, existence, and uniqueness of the
solution was done, though the numerical results exhibited some loss of contrast.
Motivated by the effectiveness of the ISS, we derive a TV framework that is globally
convex. We further construct a new ISS method for the multiplicative noise based on this
convex functional.
2. Review of multiplicative models. We are motivated by the following problem in image
restoration. In what follows we will always require f > 0, u > 0. Given a noisy image f : Ω→ R,
where Ω is a bounded open subset of R2, we want to obtain a decomposition
(2.1)
f = uη,
where u is the true signal and η the noise. We would like to denoise the signal while preserving
the maximum information about u. In image processing, such information is manifested in a
large part by keeping the sharpness of edges, since conventional denoising methods tend to
blur the image by smearing the edge information or by creating ringing artifacts at edges.
We assume that we have some prior information about the mean and variance of the
multiplicative noise,
1
N
1
η = 1,
N
(η− 1)2 = σ2,
(2.2)
where N =
to σ2.
(2.3)
1. This states that the mean of the noise is equal to 1, and the variance is equal
The procedure introduced in [22] sought the solution of a constrained optimization problem
in the space of bounded variation, implemented by the gradient projection method involving
Euler–Lagrange equations and artificial time evolution as in [23]. The following optimization
problem was studied:
u = arg min
u∈BV (Ω)
{J(u) + λH(u,f )},
Ω
|∇u| and H(u,f ) is a fidelity term for a known signal or image f . In the case
where J(u) =
of additive noise, H(u,f ) = 1
L2, whereas H(u,f ) takes a more sophisticated form for
2
multiplicative noise. A fidelity term H(u,f ) was used, consisting of two integrals with two
− 1)2. The initial data was chosen to satisfy the
Lagrange multipliers H(u,f ) = λ1
constraints
f − u2
f
u + λ2
( f
u
(2.4)
where N =
Ω 1.
1
N
1
N
f
u(0)
f
u(0)
2
− 1
= 1,
= σ2,
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 12/01/12 to 125.71.231.223. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
INVERSE SCALE SPACE FOR MULTIPLICATIVE NOISE
297
A gradient projection method was introduced to make sure the two constraints are always
satisfied during the evolution, which reduces to requiring
(2.5)
and
(2.6)
∂
∂t
=−
f
u
2
∂
∂t
− 1
f
u
= 0 =
∂
∂t
f
u2 ut = 0
f 2
u3 ut.
2
=−2
f
u
Thus we evolve the following Euler–Lagrange equation:
(2.7)
ut =∇· ∇u|∇u| + λ1
f
u2 + λ2
f 2
u3 ,
where the values of λ1 and λ2 can be found dynamically using the gradient projection method.
However, for this to be a convex model we need λ1,λ2 ≥ 0, which is not necessarily the case
as t evolves. Moreover, if we fix λ1,λ2 > 0, then the corresponding minimization problem will
lead us to a sequence of constant functions u approaching +∞.
Multiplicative noise, in various imaging modalities, is often not necessarily Gaussian noise.
In the speckle noise model, such as SAR imagery, the noise is treated as gamma noise with
mean equal to 1. The distribution of the noise η takes the form of g(η),
(2.8)
g(η) =
LL
Γ(L)
ηL−1 exp(−Lη)1{η>=0}.
Based on this, Aubert and Aujol [2] formulated the following minimization problem:
(2.9)
u = arg min
u∈S(Ω)
J(u) + λ
logu +
,
f
u
|∇u|, and new fitting function H(u,f ) =
Ω
where J(u) =
u ) is strictly convex for
u∈ (0,2f ). The derivation of this functional is based on MAP on p(u|f ), assuming the noise η
follows a gamma law with mean 1, and p(u) follows a Gibbs prior. Some interesting analysis
of this model was done in [2].
(logu + f
For denoising, the ISS method has so far been applied only to the additive noise model.
In this paper we will present a new TV functional for the multiplicative noise model, which
is globally convex. We will also extend the TV minimization approach to the ISS for our new
model.
3. Inverse scale space. In this section, we review some key results for ISS. An iterative
algorithm based on TV was introduced in [21], which was preceded by several interesting and
related pieces of work [16, 25, 26]. The ISS is based on the iterative algorithm, by taking the
time step to the limit of 0, formally introduced in [5, 6] for dealing with additive noise and
deconvolution. We review some key results for the ISS in this section.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 12/01/12 to 125.71.231.223. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
298
JIANING SHI AND STANLEY OSHER
Following the treatment of Burger et al.
[5, 6], the TV-based minimization (2.3) can
be extended into an iterative algorithm. Introduced in [21], the authors propose solving a
sequence of uk,
(3.1)
uk = arg min
u∈BV (Ω)
{D(u,uk−1) + λH(u,f )},
where uk is the primal variable and D(·,·) is the Bregman distance related to a differentiable
functional J : Rn → R, formally defined by
(3.2)
J (u,v) = J(u)− J(v)−u− v,p,
Dp
p∈ ∂J(v).
The initial condition u0 is a constant satisfying
∂uH(u0,f ) = 0 and p0 = 0. Here ∂uH(u,f )
is the subgradient of H(u,f ) with respect to u, and p∈ ∂J(u) is an element in the subgradient
of J at u. We will also use the notation pk, which is an element of the subgradient of J at uk.
For a continuously differentiable functional there exist a unique element p in the sub-
differential and, consequently, a unique Bregman distance. For a nonsmooth and not strictly
convex functional such as TV, one can still obtain convergence of the reconstruction, as long
as the functional has suitable lower semicontinuity and compactness properties in the weak-*
topology of BV (Ω) (and also in L2(Ω) by compact embedding).
One can rewrite (3.1) and arrive at
(3.3)
uk = arg min
u∈BV (Ω)
{J(u)−u,pk−1 + λH(u,f )}.
The Euler–Lagrange equation for (3.3) is
pk − pk−1 + λ∂uH(uk,f ) = 0.
Such an iterative regularization is taken to the limit by letting λ = Δt, uk = u(kΔt), and then
dividing by Δt and letting Δt→ 0. This leads to a continuous flow in the ISS,
(3.4)
∂tp =−∂uH(u,f ),
p∈ ∂J(u),
with the initial conditions p0 = 0, u0 = c0, where c0 is a constant satisfying
This was introduced in [21] and analyzed in [4].
∂uH(c0,f ) = 0.
This elegant formulation of the ISS flow is not straightforward to compute in two or more
dimensions, since the relation between p and u is complicated for TV and other nonlinear
cases. However, in one dimension one can solve this directly by regularizing the TV term
J(u) by the following:
|∇u|2 + ,
(3.5)
J(u) =
Ω
where > 0 is a small constant.
A relaxed form of the ISS flow was introduced, which involves the following coupled
equations:
(3.6)
∂tu = −p + λ(−∂uH(u,f ) + v),
∂tv = −α∂uH(u,f ),
with p∈ ∂J(u) and initial conditions u0 = c0, v0 = 0.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 12/01/12 to 125.71.231.223. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
INVERSE SCALE SPACE FOR MULTIPLICATIVE NOISE
299
Some examples of the ISS flows using different choices for J(u) and H(u,f ) are as follows:
• Linear model: J(u) = 1
∇u2
2
2, H(u,f ) = 1
2
f − u2
L2.
ut = Δu + λ(f − u + v),
vt = α(f − u).
|∇u|, H(u,f ) = 1
2
f − u2
L2.
Ω
(3.7)
• ROF model: J(u) =
(3.8)
• T V -L1 model: J(u) =
(3.9)
Ω
ut = ∇· ∇u|∇u| + λ(f − u + v),
vt = α(f − u).
|∇u|, H(u,f ) =f − uL1.
ut = ∇· ∇u|∇u| + λ(sign(f − u) + v),
vt = αsign(f − u).
Note that H(u,f ) is not strictly convex or smooth here, and sign is the notation for
an element in the subgradient of H(u,f ).
• Deconvolution by ROF: J(u) =
|∇u|, H(u,f ) = 1
f − K ∗ u2
L2.
2
Ω
(3.10)
ut = ∇· ∇u|∇u| + λ
vt = α ˆK ∗ (f − K ∗ u).
ˆK ∗ (f − K ∗ u) + v
,
Here K is a real blurring kernel, ˆK(x,y) = K(−x,−y), and ∗ denotes convolution [10].
So far (3.8)–(3.10) have not been studied analytically. However, a nice a priori inequality
was obtained in [18].
4. New methodology. The difficulty of extending previous TV models for multiplicative
noise to the ISS method is due to the lack of global convexity. We will overcome this difficulty
via the following approaches.
4.1. Logarithm transform. The simplest idea is to take the log of both sides of (2.1),
logf = logu + logη,
which essentially converts the multiplicative problem into an additive one. The additive
noise problem has already been successfully treated using the ROF method [23], the Bregman
iterative method based on the ROF model [21], as well as with the ISS and its relaxed version
[5, 6].
Now we consider seeking w = logu based on the noisy observation logf . We obtain the
following TV minimization, using the BV norm of w. We no longer use the BV norm of u,
and this gives us a convex optimization problem:
(4.1)
w = arg min
w∈BV (Ω)
J(w) +
,
w− logf2
L2
λ
2
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 12/01/12 to 125.71.231.223. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
300
JIANING SHI AND STANLEY OSHER
where J(w) =|w|BV and the fidelity term H(w,f ) = 1
L2. Note that a nice property
of H(w,f ) is its strict convexity. It follows that the corresponding relaxed inverse scale space
(RISS) flow can be expressed as
w− logf2
2
(4.2)
wt = ∇·
vt = α(logf − w),
∇w|∇w|
+ λ(logf − w + v),
with v(0) = 0, w(0) = c0, for c0 =
logf /
1.
This method gives excellent results, as we shall see in section 6. However, we will also use
a more general model, described as follows, which gives equally excellent results.
4.2. Convex total variation minimization for a general multiplicative noise model. We
consider a rather general TV formulation for the multiplicative noise:
(4.3)
u = arg min
u∈BV (Ω)
J(u) + λ
a
+ clogu
,
2
f
u
+
b
2
f
u
because a and b might be negative in [22].
where J(u) is the TV of u and a,b,c are nonnegative constants. This formulation seems to
include all previous models [2, 12, 22].
• When c = 0, this resembles the (constrained) model of [22]. It is generally different
• When b = 0 and a = c, this reduces to the model in [2].
• When a = 0 and b = c, this reduces to the model in [12].
We can solve this TV problem by gradient descent, resulting in the following unsteady
Euler–Lagrange equation:
(4.4)
ut =∇· ∇u|∇u| + λ
a
f
u2 + b
f 2
u3
− c
1
u
.
We require the following conditions:
• The fidelity term H(u,f ) has a minimum at u = f , which implies that c = a + b.
• The fidelity term H(u,f ) is convex for u near f , which means a + 2b > 0. However,
this term will not be convex for u f .
Due to the lack of global convexity, we abandon this approach and proceed by letting
w = log(u). We also replace J(u) = J(exp(w)) by J(w). By using such a partial transformation
of variables, the fidelity term H(w,f ) is rendered convex, while the TV term is left unaltered.
Such replacement is reasonable since the mapping w→ ew is monotically increasing. Convexity
of the objective function allows us to take advantage of the excellent regularization ability of
the inverse scale space. It is noteworthy that we use the BV norm of w, instead of u. The
regularization of the BV norm therefore takes place in the logarithmic sense for the image,
penalizing on log(u) instead of u.
In this logarithmic case, the amount of regularization
depends on the image intensity, resulting in strong smoothing for image values near 0, and
less smoothing for large intensities. This is somewhat visible in Figures 1 and 2.
One can obtain the following TV based minimization:
(4.5)
w = arg min
w∈BV (Ω)
J(w) + λ
af exp(−w) +
b
2
f 2 exp(−2w) + (a + b)w
.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 12/01/12 to 125.71.231.223. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
INVERSE SCALE SPACE FOR MULTIPLICATIVE NOISE
301
The fidelity term H(w,f ) =
We aim to solve this minimization problem for w and recover the signal by u = exp(w).
convex. Using gradient descent and the Euler–Lagrange equation for this TV based problem
we obtain
2 f 2 exp(−2w) + (a + b)w) is globally strictly
(af exp(−w) + b
af exp(−w) + bf 2 exp(−2w)− (a + b)
.
(4.6)
wt =∇· ∇w|∇w| + λ
We choose the initial condition w(0) defined below so that
∂wH(c0,f ) = 0.
In order to perform iterative regularization on this model, we define the following sequence
{wk} together with {pk}.
• Set w0 = c0, p0 = 0.
• Given wk−1 and pk−1 ∈ ∂J(wk−1), k≥ 1, perform the following two steps:
(i) Compute wk = argmin
Qk(w) with
w
Qk : w−→ J(w)− J(wk−1)−pk−1,w− wk−1 + λH(w,f ),
(af exp(−w) + b
2 f 2 exp(−2w) + (a + b)w) and ·,· denotes the
where H(w,f ) =
usual duality product.
(ii) Update the dual variable
pk = pk−1 + λ∂wH(w,f )∈ ∂J(wk),
where ∂wH(w,f ) =−af exp(−w)− bf 2 exp(−2w) + (a + b).
• Increase k by 1 and continue.
Such an iterative procedure can be further extended to the nonlinear ISS flow; see (3.4).
We obtain the following ISS model:
= af exp(−w) + bf 2 exp(−2w)− (a + b),
∂p
∂t
p ∈ ∂J(w),
For the initialization of the ISS, w(0) = c0, we use the fact that
(4.7)
with w(0) = c0, p(0) = 0.
(4.8)
∂tp = 0,
Ω
which means that, for all time,
t→ 0, leading to
Ω ∂wH(w(t),f ) is time invariant. This expression vanishes as
− af exp(−c0)− bf 2 exp(−2c0) + (a + b)
⎧⎪⎪⎪⎪⎪⎨
⎪⎪⎪⎪⎪⎩
1
⎛
⎝−a
−log
−log
a + b
(a
f )2 + 4b
f 2
f +
f
2b
f 2
= 0.
if b = 0,
⎞
⎠ if b= 0.
This gives us the initial condition for the ISS flow:
(4.9)
c0 =
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 12/01/12 to 125.71.231.223. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
