Reconstruction from 2-D wavelet transform
modulus maxima using projection
A.W.C.Liew and N.F.Law
Abstract: Wavelet transform modulus maxima can be used to characterise sharp variations such
as edges and contours in an image. The authors analyse the a priori constraints present in the
wavelet transform modulus maxima representation. A new projection-based algorithm which
enforces all the a priori constraints in the representation is proposed. Quadratic programming is
used to obtain a sequence which satisfies the maxima constraint, thus realising the projection onto
the maxima constraint space. To save computation, an approximate method to obtain a sequence
which satisfies the maxima constraint is given. The new algorithm is shown to provide better
solution than the original reconstruction algorithm of Mallat and Zhong. The authors also propose
a simple method to accelerate the algorithm. The acceleration is achieved by the incorporation of a
momentum term which exploits the high correlation between the difference images between two
consecutive iterations. The simulation results show that the proposed algorithm gives good
reconstruction and the simple acceleration method can significantly improve the convergence rate.
1 Introduction
Points of sharp variation such as edges and discontinuities
are usually one of the most important features for analysing
the properties of signals or images. However, since these
features usually exist in a range of scale, a multiscale
strategy needs to be adopted. The need to adopt a multi-
scale strategy for describing a signal or image is also
brought out by the study on the biological visual system.
It was conjectured that the basic representation, the primal
raw sketch, furnished by the retinal system is a succession
of contour sketches at scales that are in geometric progres-
sion [1]. The wavelet transform modulus maxima repre-
sentation (WTMM) proposed by Mallat and Zhong [2]
provides such a multiscale contour representation of an
image. This representation is obtained by retaining the
local modulus maxima of the continuous dyadic wavelet
transform (WT) which correspond to discontinuities in an
image. It provides a compact but physically meaningful
description of an image. Two important issues related to
the WTMM are: first, 1. the uniqueness of the representa-
tion; and secondly, 2. the reconstruction of an image from
its WTMM. The uniqueness issue of the representation has
been investigated for 1-D signal [3, 4]. The uniqueness of
the representation depends on the total number of WT
modulus maxima in the signal. If the maxima points are
sparse, the representation is, in general, not unique. In that
IEE, 2000
IEE Proceedings online no. 20000206
DOI: 10.1049/ip-vis:20000206
Paper first received 19th October 1998 and in revised form 30th November
1999
A.W.C. Liew is with the Department of Electronic Engineering, City
University of Hong Kong, Tat Chee Avenue, Kowloon Tong, Hong Kong
N.F. Law is with the Department of Electronic and Information
Engineering, Hong Kong Polytechnic University, Hung Hom, Kowloon,
Hong Kong
176
situation, it will still be desirable if a close approximation
of the original signal can be reconstructed.
A number of reconstruction methods have been
proposed for the 1-D WTMM reconstruction problem
such as the projection-based method [2, 5], the conjugate
gradient error minimisation method [6] and the least square
eigenspace method [7]. For the 2-D WTMM reconstruction
problem, Mallat and Zhong [2, 8] proposed a projection
onto convex space method. In this method, the a priori
constraints are formulated as some convex sets. An alter-
nating projection technique is then used to obtain a feasible
solution at the intersection of all the convex sets. As one of
the a priori constraint in the WTMM is not convex, their
algorithm does not enforce the constraint strictly [4]. We
will show that this affects the quality of the reconstruction.
In this paper, we propose a new projection algorithm
which enforces the a priori constraint strictly in an efficient
way. For any iterative algorithm, a main concern is the
convergence rate. We propose a way to improve the
convergence rate based on the idea of momentum. Some
experimental results will be presented to show the perfor-
mance of our algorithm.
2 The reconstruction problem
Let the wavelet functions (x, y) and Y²(x, y) be the
partial derivatives of a 2-D smoothing function, i.e. a cubic
spline [2,8]. Thus the 2-D dyadic WT of an image ƒ(x, y)
becomes the gradient of the image at multiples of dyadic
scales. The two components of the 2-D dyadic WT are
denoted as W, f(x, y) and W2,f(x, y). The local modulus
maxima of W½, f(x, y) along x, for constant y, are the points
of sharper horizontal variation of f(x, y) smoothed at the
scale 2. Similarly, the local modulus maxima of W2, f(x, y)
along y, for constant x, are the points of sharper vertical
variation. The local maxima belong to curves in the (x, y)
plane which are the edges of the image along each
IEE Proc.-Vis. Image Signal Process., Vol. 147, No. 2, April 2000
direction [2, 8]. Let X be the set of modulus maxima
locations of W};f(x, y) along the horizontal direction, i.e.
X = { xk+ Yk \ Wy f(x, y) has local maxima at
(1)
(x, y) = (x, y) along the horizontal direction}
The collection of X for all scales j is denoted by
X₁ = {X}jez. Similarly
X} = {x1,y1 : \W}; ƒ (x, y)] has local maxima at
(x, y) = (x,y) along the vertical direction}
(2)
and X2 = {X}};cz. For a J-level 2-D dyadic WT, the
collection
{W}ƒ (xk+ Yk); W}ƒ (x1, y₁), S21ƒ (x, y)}
(3)
where 1≤j≤J, (Xk, Yk) ¤ X 1, (X1, y₁) € X2, is called the 2-D
WTMM of f(x, y). S₂f(x, y) is the lowpass approximation
at level J.
For the 2-D WTMM of an image ƒ to be unique, the set
of wavelets {(Xks Yk), P2/(x1, 11)} that corresponds to
the maxima locations has to constitute a frame in the 2-D
signal subspace, i.e. there exists two constants A, B>0,
A<B< ∞ such that
A\\ ƒ ² ≤ Σ Σ \ < ƒ, ¥½ (xx. 3x) > β [ !
jez \(xkk)
\ <ƒ, y
(x1,91)
< ƒ, ¥2; (x1,y1) > |² ) ≤ B||ƒ||² (4)
If eqn. 4 is satisfied, then an image can be reconstructed
exactly from its WTMM. While it is possible to determine
the uniqueness of a 1-D WTMM [3, 4], it is unclear as to
how to determine the uniqueness of a 2-D WTMM in
practice. A more practical issue, however, is related to the
problem of how to reconstruct an image from its WTMM.
The 2-D reconstruction attempts to reconstruct the
image by first reconstructing its 2-D dyadic WT from its
2-D WTMM. The image is then obtained by performing a
2-D inverse dyadic WT on the recovered transform.
Let K be the solution space of the 2-D reconstruction
problem, then any element (g(x, y), g(x, y)}jez belonging
to K must satisfy the following constraints:
1. It must be the 2-D dyadic WT of a function in L2(R²)
2. At each scale 2/ and for each modulus maxima location
(x, y) EX and (x, y) EX2, we have
27
g}(xk+ Yk) = W₂f (xk+ Yk)
It is clear that imposes constraint 1 and I imposes
constraint 2. The maxima constraint 3 is, however, not
imposed by either For V. There is no guarantee that an
element in A will only have modulus maxima at the
prescribed locations. Any signal having its dyadic WT
passing through those wavelet modulus maxima points is
by definition in A and is considered a mathematically valid
solution. One can easily see that the true solution space K
is a subspace of A.
To 'approximately' impose constraint (iii) into the solu-
tion, Mallat and Zhong required that: the energy of the
reconstructed signal's dyadic WT be as small as possible
on average (which will tend to create local modulus
maxima only at the prescribed locations) and that the
energy of the derivative of the reconstructed signal's
dyadic WT be as small as possible on average (which
will tend to suppress local modulus maxima at other
locations). The above considerations led to the proposal
of minimising the Sobolev norm when projecting onto the
space г. The projection operator Pr is defined to be the
operator that transforms any sequence of functions {g}(x,
y), g}(x, y)} jez onto the sequence of functions {h}(x, y),
h(x, y)}jez ET closest to it with respect to the Sobolev
norm. Mathematically, this is expressed as
Pr : (g}, 8});ez → (h}, h} )jez
wwww..com
(7)
such that the Sobolev norm of (}, })jez=(h} — g},
h} - g})jez, i.e.
∞
11½ 11² + 1 &² 11² + 2²²
(
θε
+
Эх
ду
(8)
is minimised. Their algorithm is then stated as iterating the
composite operator PPP, i.e. y lim∞P"x, where x
and y are the initial estimate and the final solution,
respectively. Although Pr will tend to create modulus
maxima at the prescribed locations while trying to suppress
spurious maxima at other locations, the projection onto г
is no longer a minimum distance projection with respect to
the Euclidean norm L₂ and therefore not necessarily a
nonexpansive operator with respect to the L₂ norm. The
iterative algorithm therefore has the potential to diverge in
L₂ norm. In addition, constraint 3 is not necessarily
satisfied upon convergence.
Although we found that the reconstruction algorithm of
Mallat and Zhong generally produces satisfactory solution
in most cases, poor reconstruction can sometime arise. To
illustrate, we show in Figs. 1-4 1-D reconstruction of a
signal from its WTMM using Mallat and Zhong's
g}(x1,y1) = W}ƒ (x1, Y1)
(5)
3. At each scale 2/, the 2-D WT modulus maxima obtained
from g(x, y) and g(x, y) are located at the abscissa
(X, Y)EX and at the abscissa (x, y) EX2, respectively.
The 2-D projection based reconstruction algorithm [2, 8]
employs alternating projections onto an affine space à
specified by the modulus maxima and onto the 2-D
dyadic wavelet space V The space I is defined as the
affine space of the sequence of functions {g}(x, y), g}(x,
y)}jez such that for all scales 2 and for all modulus
maxima locations (xk, уk) €X] and (x1, y₁) €X½ eqn. 5 is
satisfied. The 2-D dyadic wavelet space V is defined as the
space of all sequences of functions that satisfy the 2-D
reproducing equation. The solution to their reconstruction
thus lies on the intersection of ½ and Ã, i.e.
1.0
0.8
0.6
0.4
A: For
IEE Proc.-Vis. Image Signal Process., Vol. 147, No. 2, April 2000
0.2
0
-0.2
0
20
40
60
80
100
120
iterations
Fig. 1
original signal
Reconstruction using P
reconstructed signal obtained after 80 iterations where its SNR is at a
(6)
maximum
177
A
1
O
LM
-1
-2
2
-
W2
O
-1
-2
SNR, dB
25
20
w3
w4
54
2
1
0
-1
-2
2
1
O
-1
-2
0.8
0.4
0
-0.4
20
40
60
iterations
80
100
120
Fig. 2 Discrete dyadic WT and WTMM (impulses in wl to w4 and the
entire s4) of the original signal in Fig. 1
LM
w2
W3
W4
$4
A
-1
2
2
1
-1
-1
-2
O
+ 2
0.8
0.4
0
-0.4
0
20
20
40
60
iterations
80
100
120
Fig. Discrete dyadic WT of the reconstructed signal in Fig. 1
178
15
0
20
40
60
iterations
80
100
120
Fig. 4 SNR against number of iterations using P for the reconstructed
signal of Fig. 1
algorithm. For this signal, its WTMM is not unique and
exact reconstruction is not possible. The iterative recon-
struction is started from the zero sequence. The recon-
structed signal shown in Fig. 1 is of a fairly poor quality
especially for the square pulse. Fig. 4 shows that Mallat
and Zhong's algorithm exhibits oscillatory behaviour and
fails to converge monotonically in L₂ norm. By comparing
Fig. 2 with Fig. 3, we observe that spurious maxima are
present in the WTMM of the final reconstructed signal.
The dyadic WT of the reconstructed signal oscillates
slightly at locations of abrupt change. This is in agreement
with the observation reported by Mallat and Zhong [2], in
which they proposed a heuristic way to clip these spurious
maxima.
3 The proposed solution
We propose a projection-based reconstruction algorithm
that incorporates all the available a priori constraints into
the reconstruction such that the solution must only have
WT modulus maxima at the locations X₁ and X2.
The projection operator onto the linear 2-D dyadic
wavelet space that imposes constraint 1 is given by the
reproducing kernel equation
Py = WW-1
(9)
Let C be the space that imposes constraints 2 and 3, i.e.
Cr Space for constraint 3. Then the projection Pc can
be implemented as follows. First, the values of the esti-
mated 2-D dyadic WT at locations X, and X₂ are replaced
by the true values. We could then have two possible
situations,
1
1 A signal which is not monotonically increasing or
decreasing between two modulus maxima of opposite
sign. A monotonically increasing or decreasing signal
that is closest (in L₂ sense) to the given signal needs to
be constructed.
2 A signal which has more than one local modulus mini-
mum between two adjacent modulus maxima of the same
sign. A closest signal that contains only one modulus
minimum needs to be constructed. This is done by break-
ing the original signal into two signals at the point with
smallest magnitude and processing the two signals into a
non-increasing signal and a non-decreasing signal.
Notice that in situation 2 above, the point at which to break
the signal changes from iteration to iteration. This strategy
avoids the need to pre-determine a local modulus minimum
point, which is not known a priori. Also, it ensures that the
new signal is the closest to the original signal in L₂ norm
IEE Proc.-Vis. Image Signal Process., Vol. 147, No. 2, April 2000
while enforcing the modulus maxima constraint at each
iteration.
Since a separable 2-D dyadic WT is used, the 2-D
implementation of Pc is simplified to two 1-D implemen-
tations; one on the rows of Wg(x, y) and the other on the
columns of W2g(x, y). In discrete implementation, the
projection Pe can be formulated as a quadratic program-
ming problem as shown next.
N1
3.1 An exact implementation of Pc by QP
Suppose we have a discrete sequence {k}i=1,.
between the positive modulus maximum k。 and the nega-
tive modulus maximum ky that violates the nonincreasing
constraint, then the minimum distance projection of
{k;} i=0,.
i=0,………î onto C is equivalent to finding a nonincreas-
ing sequence {x;}i=0,………,n with xoko and xy=ky that is
closest to {k};=0, N · in l½ norm. The l½ distance measure
between the sequence {x;} and {k} is given by
N
The KT conditions transform the QP into a set of linear
equations involving equalities
Cj − 2 Σ 9 ;kxx - Σ ¿¡¤¡¡ +µ; = 0 \j
k
Σajx;+s;= b; ¥i
(14)
where λ;s;=µ;; =0, λ¡, µj, Xj, §; ≥0Vi, j which can be
solved using the simplex method [9]. The composite
operator for our algorithm is then given by P₁ =PyPc·
Figs 5-7 shows the improvement in the 1-D reconstruction
using our algorithm. Convergence is evident and there is
no spurious maxima in the WT of the reconstructed signal.
1.0-
0.8
N-1
d²
=
Σ(x; - k;)²
(10)
0.6
0.4
Minimisation of d² is equivalent to maximisation of
N-1
i=1
N-I
- ΣNTx+2Σx. With no loss of generality, we
can add a constant term to the sequence {k};-0, ..
that k₁ > 0, k = 0. Then we have the following discrete
time optimisation problem:
Given ko, k₁,...,ky, find x¡,
0.2
N SO
0
-0.21
0
20
40
N - 1 to maxi-
60
iterations
80
100
120
mise
N-1
N-1
x² + 2 Σ kpx;
(11)
i=1
subject to
Xj ≤ ko
-
.., N
-
2
X;+1 − x; ≤ 0_ i=1,...,
x; ≥ 0
i = 1,..., N − 1
Since the objective function is quadratic and positive
semidefinite and the constraints are linear, it is a concave
quadratic programming (QP) problem [9]. For a concave
QP problem, the optimal solution always exists and satis-
fies the Kuhn-Tucker (KT) conditions:
Given a QP
max
ΣΣ
Cjxj − Σ xjQjk Xk
j,k
subject to apx; ≤b; Vi
Σ
Vi
(12)
W4
W3
Fig. 5 Reconstruction using P
original signal
reconstructed signal obtained after 100 iterations
W2
w1
-1
-2
2
1
1
O
O
-1
-2
2
-
M
A
-1
-2
2
1
0
-1
~
KT conditions
c; - 2Σ9;kî kΣ ha₁ + µ; = 0 vj
Vi, μ≥0
Vj
%; ≥0
λ;S; =0 Vi, μ;x=0 Vj
where {x;}, {$;} are the optimal solution to QP.
IEE Proc.-Vis. Image Signal Process., Vol. 147, No. 2, April 2000
0.8
0.4
0
(13)
-0.4
0
20
40
60
80
100
120
iterations
Fig. 6 Discrete dyadic WT of reconstructed signal in Fig. 5
179
SNR, dB
40
30
20
1
20
40
60
80
100
iterations
*
*
LA
*
Fig. 7 SNRs against number of iterations using P (solid line) and P¡
(dotted line)
3.2 An approximate implementation of Pc
Unfortunately, the quadratic programming technique is
computationally too expensive in a 2-D setting. Instead,
we describe a simple approximation that implements the
projection Pc.
Consider the signal shown in Fig. 8. This signal has one
spurious maximum in between two WT modulus maxima
m₁ and m2. The problem is to find a non-increasing signal
between m₁ and m₂ that is closest to the given signal in L2
norm. Obviously, the new signal will be exactly the same
as the given signal except at the interval between x2 and x5.
At this interval, the signal should be replaced by a constant
value v such that the sum of the two L₂ errors e₁ and e2 is
minimised. We approximate the above problem to finding
an optimum combined estimate v from v, and v₁ that has a
minimum weighted error. Let E and E₂ be the L₂ error
associated with choosing v=v₂ and v=v, respectively.
The two errors E, and E2 are given by
E\
=
(f(x) – v₁)² dx
m₁
D
*
m2
Fig. 9 Obtaining a nonincreasing signal between two modulus maxima
mĮ and mz when spurious maxima are present
ar
E2
= [ "
(f(x) — v₁)² dx
(15)
xz
If E₁ is much smaller than £2, then v should be much closer
to v than to v i.e. the term (v - v₂)² should be much
smaller than (vv)². This would be achieved if, in
finding the least sum of the two terms, we weigh the
former term by a larger constant than the latter term. A
f(x)
V1
VD
V
V
a
V2
given signal
new signal
컴터
NZ NZ
E2
மிக
02
M2
m1
X1 X2
X3
X4 X5X6
X
Fig. 8 Signal with spurious maximum between modulus maxima
m2
180
Fig. 10 Original ‘Lena' image, 256 × 256 pixels
m¡
and
Fig. 11 Original 'Peppers' image, 256 × 256 pixels
IEE Proc.-Vis. Image Signal Process., Vol. 147, No. 2, April 2000
Fig. 12 Four-level dyadic WT and WTMM of 'Lena'
First two columns from the left show, respectively, Wƒ and W3, f for scale 1 to 4
Third and fourth columns show the WT modulus maxima obtained from W,falong the horizontal direction and from W3 ƒ along the vertical direction, respectively
Bottom image is the lowpass approximation at level 4
logical weighting is to weight each term by one-over their
associated error
-
E
(v − v₂)² __ (v − v¿)²
E₂
E₁
+
(16)
In eqn.
16, the two terms are automatically weighted
according to their importance. Differentiating eqn. 16
with respect to v and setting the result equal to zero yields
y = Va
E₂
+
-(1%)
E₁ + E₂
Va)
(17)
When a signal has more than one spurious maximum, it
can be processed as illustrated by the example in Fig. 9 to
give a signal with no spurious maxima. A spurious maxi-
mum can only be removed if the interval at which E₁ or E₂
is calculated contains only one minimum or one maximum.
This ensures that only one spurious maximum is removed
at a time. The magnitude of a signal between the interval
m₁ and
m2 is also required to lie between V₁ and
V2
before
any spurious maximum is smoothed out. This can be done
by replacing any value that is larger than v₁ by v₁ and any
value that is smaller than v₂ by v2. The following steps are
taken when removing multiple spurious maxima:
IEE Proc.-Vis. Image Signal Process., Vol. 147, No. 2, April 2000
181
1. Remove the left most spurious maximum (marked by an
asterisk in Fig. 9) first if possible
2. If the left most spurious maximum cannot be removed,
move to the next spurious maximum on the right until a
spurious maximum is removed
3. Repeat 1 and 2 until all spurious maxima are removed
In actual implementation, the above procedure is
performed on discrete time samples. The L₂ norm becomes
the norm and the integral becomes summation. The
discrete sequence obtained from the above procedure
becomes an approximation to the optimum discrete
sequence obtained using QP. The approximate sequence
was observed to provide a very good approximation to the
optimum sequence obtained using QP. Simulation with 1-D
signals indicates that reconstruction using the approximate
sequence is almost identical to that using the exact imple-
mentation of Pc via QP. While our projection algorithm
using the exact implementation of Pe is computationally
more expensive and hence slower than Mallat and Zhong's
projection algorithm, our projection algorithm using the
approximate Pc is comparable in computation speed to
their projection algorithm.
a
a
b
Reconstruction 'Lena' from four-level subsampled WTMM after
Fig. 13
50 iterations
a Reconstruction using P
b Reconstruction using P
40
10
SNR, dB
30
35
36
25
25
20
20
0
20
40
60
80
100
iterations
b
Fig. 15 Reconstruction of 'Peppers' from four-level subsampled WTMM
after 100 iterations
a Reconstruction using P
b Reconstruction using P,
SNR, dB
45
40
55
35
30
25
20
20
0
20
40
60
80
100
iterations
Fig. 14 Reconstruction of 'Lena' from four-level subsampled WTMM
after 50 iterations
Fig. 16 Reconstruction of ‘Peppers' from four-level subsampled WTMM
after 100 iterations
SNR against number of iterations
p
Pi
182
SNR against number of iterations
Р
P
IEE Proc.-Vis Image Signal Process., Vol. 147. No. 2, April 2000
4 Simulation results
To compare our proposed algorithm with Mallat and
Zhong's algorithm, we present some reconstruction results
for two images shown in Figs 10 and 11. Fig 12 shows the
four-level dyadic WT and WTMM of 'Lena'. Since the
lowpass approximation has most of its energy concentrated
in a bandwidth that is approximately 16 times smaller than
the bandwidth of the original signal in the horizontal and
vertical directions, it can be subsampled by a factor of 24 in
both directions, i.e., the 256 × 256 image is subsampled to
16 × 16. To reconstruct the image from its subsampled
WTMM, the subsampled lowpass approximation needs to
be upsampled back to the original dimension first. In these
examples, simple linear interpolation between samples is
used for upsampling.
The reconstructed images of 'Lena' obtained after 50
iterations and the reconstructed images of 'Peppers'
obtained after 100 iterations using P and P₁ are shown in
Figs 13-16. For the image 'Lena', both P and P₁ gave
results that are visually identical to the original image after
50 iterations. The divergence of Mallat and Zhong's algo-
rithm in L₂ norm is evidence in the reconstruction for
'Lena'. For the image 'Peppers', the reconstruction using
P₁ has better contrast than the reconstruction using P. The
reconstruction using P has better convergence initially and
this is because the operator Pr in P projects to a sequence
that is closer to the true sequence initially.
Γ
5 The accelerated algorithm
It is possible to accelerate the convergence of the projec-
tion algorithm P₁ by exploiting the similarity of the
1
difference image between each iteration. Let f be the
reconstructed image at the ith iteration, the difference
image d is defined by
d; = fi+1-fi
(18)
Fig. 17 shows the difference images d1, d2, d4, ds for the
reconstruction from the subsampled WTMM of 'Lena'
using P. The difference images d and d are very similar
for both 'Lena' and 'Peppers'. A quantitative measure of
the similarity of the difference images can be obtained by
using the following similarity metrics S:
1
S(d;, di-1) =
< dj, di-1 >
||d;||||d;_1||
(19)
The value of S ranges between 0 and 1. When $= 1, d; and
d are identical to each other. Table 1 lists the value of S
between the difference images for the reconstruction from
the subsampled WTMM of 'Lena' and 'Peppers' using P₁.
The similarity between difference images d¡ and d¡— is
greater than 0.9 for i greater than 2. In fact, S approaches 1
as i increases. When S≈ 1, ƒ~ƒ¡-1+fi−1 − fi−2), i.e. fi
can be obtained from the weighted sum of f¡_¡ and ƒ¡-2
instead of from the projection P₁ on fi-1. Therefore, the
reconstruction can be accelerated by the incorporation of a
momentum term x, i.e.
1
f = f; + af; − fi-1)
1
(20)
where f and f... are the current reconstruction and the
previous reconstruction, respectively, and f'; P₁fi-1. The
momentum term a lies between 0 and 1 and it is usually
chosen to be roughly equal to S. This acceleration method
tries to build inertia into the estimation, and thus the next
estimate is predicted using the information from the
с
2
b
Fig. 17 Difference images for reconstruction of ‘Lena' from subsampled WTMM using P¡
a d₁
b dz
c da
d ds
IEE Proc.-Vis. Image Signal Process.. Vol. 147, No. 2. April 2000
Վ
183
Table 1: Similarity values of the difference images for
the reconstruction from the subsampled WTMM of
'Lena' and 'Peppers' using P₁
'Lena'
'Peppers'
S(d₂, d₁)
0.88016
0.88440
S(dz, d₂)
0.97057
0.97464
S(d4, d3)
0.98864
0.99070
S(ds, da)
0.99366
0.99480
S(d6,
ds)
0.99575
0.99636
S(d7, do)
0.99673
0.99709
S(dg, dy)
0.99727
0.99757
S(dg, do)
0.99779
0.99790
S(dio, dg)
0.99816
0.99814
a = 1 will give near optimum reconstruction performance
for most images. Fig. 18 shows the SNRs of the recon-
structions from the subsampled WTMM of 'Lena' and
'Peppers', respectively, using P₁ with a 1. The accelera-
tion effect can be readily observed.
1
The reconstruction algorithm using P₁ is relatively stable
with respect to perturbation in the WTMM. In one experi-
ment, we linearly quantised the modulus maxima and the
lowpass approximation in the WTMM of 'Lena' and
'Peppers' to 4 bits/pixel. Due to quantisation error, the
quantised subsampled WTMM of 'Lena' and 'Peppers'
have a SNR of 18.5 dB and 18.7 dB, respectively. The
reconstructed 'Lena' and 'Peppers' obtained after 20 itera-
tions have a SNR of 33.41 dB and a SNR of 31.16 dB,
respectively. The reconstructions are very accurate even
though the representations have a fairly low SNR.
previous estimate. Significant improvement in convergence
rate based on the momentum idea has also been demon-
strated [10] for blind deconvolution and in neural network
weight update [11].
For the reconstruction of ‘Lena' and 'Peppers' from their
subsampled WTMM using P₁, a > 0 gives improvement on
the unaccelerated reconstruction and ≈ = 1 gives the best
result. This is in accordance with the similarity values
presented in Table 1. In general, we found that a value of
40
35
SNR, dB
35
30
6 Conclusions
We have proposed a new projection-based algorithm to
reconstruct an image from its WTMM. In contrast to
Mallat and Zhong's algorithm, our algorithm incorporates
all the a priori constraints in the WTMM into the recon-
struction. Since one of the projection operator. Pc, which
projects a sequence into another sequence that satisfies the
modulus maxima constraint, is computationally expensive
to implement exactly using the method of QP, we provide
an approximate implementation which gives an almost
identical result. We also proposed a way to accelerate the
reconstruction algorithm based on the idea of momentum.
By exploiting the similarity of the difference images
between two consecutive iterations, we are able to use
information from the previous estimates to predict the new
estimate, thus improving the convergence rate. Simulation
results indicated that good solutions can be obtained from
our algorithm.
25
SNR, dB
20
7
20
}
20
40
60
80
100
1
iterations
A
50
15
45
40
55
35
309
30
25
20
20
O
2
References
MARR, D.: 'Representing visual information', in HANSON, A. and
RISEMAN, E.M. (Eds.): 'Computer vision systems' (Academic Press,
New York, 1978) pp. 61–80
MALLAT, S., and ZHONG, S.: 'Characterization of signals from
multiscale edges', IEEE Trans. Pattern Anal. Mach. Intell., 1992, 14,
(7), pp. 710-732
3 LIEW, A.WC., and NGUYEN, D.T.: 'Uniqueness issue of the wavelet
transform modulus maxima representation and a least square reconstruc-
tion algorithm', Electron. Lett., 1995, 31, (20), pp. 1735-1736
LIEW, A.W.C.: 'Multiscale wavelet analysis of edges: issues of unique-
ness and reconstruction'. PhD thesis, University of Tasmania, Australia,
1996
4
5
LIEW, A.W.C., and NGUYEN, D.T.: 'Reconstruction from wavelet
transform modulus maxima using nonexpansive projections', Electron.
Lett., 1995, 31, (13), pp. 1038-1039
6 LAW, N.F., and LIEW, A.W.C.: 'Wavelet maxima reconstruction by
conjugate gradient', Electron. Lett., 1997, 33, (23), pp. 1928-1929
7 LIEW, A.W.C., LAW, N.F., and NGUYEN, D.T.: 'A direct reconstruction
method for the wavelet transform extrema representation', IEE Proc.,
Vis. Image Signal Process., 1997, 144, (4), pp. 193–198
MALLAT, S., and HWANG, W.L.: 'Singularity detection and processing
with wavelets', IEEE Trans. Inf. Theory, 1992, 38, (2), pp. 617-643
9 HARTLEY, R.: 'Linear and nonlinear programming: an introduction to
linear methods in mathematical programming' (Wiley, Chichester, 1985)
10 LAW, N.F., and NGUYEN, D.T.: 'Improved convergence of projection
based blind deconvolution', Electron. Lett., 1995, 31, (20), pp. 1732-
1733
20
40
60
80
100
8
iterations
b
Fig. 18 SNRs of reconstructions from subsampled WTMM using P,
unaccelerated
accelerated (2: D
a 'Lena'
b 'Peppers'
wwwww
11 JACOBS, R.A.: ‘Increased rates of convergence through learning rate
adaptation', Neural Netw., 1988, 1, pp. 295-307
184
IEE Proc.-Vis. Image Signal Process.. Vol. 147, No. 2, April 2000
