Home

Characterization of Signals from Multiscale Edges - CMAP

image

Contents

1. __ _ _ i a CCC 8 9 c Fig 2 a Signal of 256 samples b discrete dyadic wavelet transform of signal a computed on nine scales At each scale 2 we plot the signal W f which also has 256 samples c modulus maxima of the dyadic wavelet transform shown in b Each Dirac indicates the position and amplitude of a modulus maximum J AP B D DAN e One can derive from this equation that the higher frequencies of S f x which have disappeared in So f x can be re covered from the dyadic wavelet transform Wz f y lt j lt J between the scales 2 and 27 We suppose that the original signal is a discrete sequence D dn nez Of finite energy If there exists two constants C gt 0 and C gt 0 such that o w satisfies 00 Ww e R C lt 5 l lw 2nr lt C2 23 n 00 714 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE VOL 14 NO 7 JULY 1992 then one can prove prove 18 that there exists a function f x L R not unique such that Yne Z Sif n dy 24 The input signal can thus be rewritten D Si f n lt 7 For a particular class of wavelets defined in Appendix A the discrete signal D Si f n lt 7 allows us to compute a uniform sampling of the wavelet transform of f a at any scale larger than 1 Let us denote Wy f Wo f n w neg and S f Sai f n w nez 25 where w is a sampling shift that depends only on y 2 For any c
2. 2 X PpX de x de xi 2 oio e a Fe x 4 i ae n j EZ 121 We know that the function x is the sum of two exponentials given by 113 between any two consecutive modulus maxima located at x2 and z 1 If we replace the constants a and 8 with their values specified by z3 and zi 1 With a few algebraic manipulations we derive that 2e x seine 275 zi pema a gt In dx Th OMO a raa 122 9j TUE e x The derivative at zf 4 is the left derivative whereas the derivative at n is the right derivative Since we suppose that there exists a constant D gt 0 such that 271 x x _ gt D lt 1 we obtain X PpX gt gt E lueh 123 which proves 119 for Cy 2 Let us now prove that there exists a constant C2 such that for any X Wzig z ez V IX Pax so SO Alae n j EZ Let U and O be the spaces defined in Section V B The function g x can be decomposed into g x gi x g2 2 with gi z U and go x O The original function f s can also be decomposed into f x file fala with f x U and fo x O Let us now define the function h x fi x g2 2 127 Since h x f x u z with u x O we know from 46 that h x satisfies constraint 1 and thus that Wash I We also have h x g x fi x gi x U Since we suppose 124 125 126 732 IEEE TRANSACTIONS ON PATTERN ANALYSIS
3. K Zwe 2wy e K Wy L we b We Wy where the functions K w and L w are 27 periodic and satisfy G w K w H w 1 107 1 H w aaa Gee One can prove that wz wy and X7 w wy are the Fourier transform of reconstructing wavelets that satisfy 61 As in Appendix A we choose G w 4ie sin w 2 to approximate a derivative The 2 D wavelets used in the computations of this article are derived from the 1 D quadratic spline wavelet shown in Fig 1 The values of the discrete filters H G K and L are given in Table I L w 108 APPENDIX D FAST WAVELET ALGORITHMS FOR 2 D SIGNALS We describe two fast algorithms to implement the wavelet transform and the inverse wavelet transform in two dimen sions We suppose that the two wavelets y x y and Yy x y are characterized by the three discrete filters H G K and L mentioned in Appendix C We use the same notations as in Appendix B and L is the discrete filter obtained by putting 2P 1 zeros between consecutive coefficients of the filter L We also denote by D the Dirac filter whose impulse response is equal to 1 at 0 and 0 otherwise We denote by A H L the separable convolution of the rows and columns respectively of the image A with the 1 D filters H and L 730 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE VOL 14 NO 7 JULY 1992 The following algorithm computes the 2 D discrete wavelet transform of an image S f A
4. at the scale 27 has two components defined by Wa f a y f Yz x y and W3 f x y f Yz 2 y 56 MALLAT AND ZHONG CHARACTERIZATION OF SIGNALS 60 50 40 30 20 10 0 200 400 600 800 1000 Fig 7 SNR of the reconstructed signal computed with respect to the signal to which we converge as a function of the number of iterations on the operator P We refer to the 2 D dyadic wavelet transform of f x y as the set of functions Wi Wa f x y W f 2 9 jez Let We Wy and b We Wy be the Fourier transforms of p x y and Y x y The Fourier transforms of W4 f x y and W3 f x y are respectively given by 57 Wy J We Wy F we wyj Z ws 2 wy 58 We f We Wy flws wyj Z wr 2 wy 59 To ensure that a dyadic wavelet transform is a complete and stable representation of f x y we impose that the 2 D Fourier plane is covered by the dyadic dilations of y1 w wy and Wp we Wy This means that there exist two strictly positive constants A and B4 such that V We Wy R Lay aa ale a Ag lt rcs we wy J 60 Let x x y and x x y be two functions whose Fourier transform satisfy co 2 j Hee wy 2 we 2 wy GP Were Hisy we Py 1 61 There is an infinite number of choices for x x y and x x y We can derive from 58 59 and 61 that f x y is reconstructed from its dyadic wavelet transform with oo f x y W
5. then O 0 which implies that h x must be equal 46 718 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE VOL 14 NO 7 JULY 1992 to f z In general this is not the case therefore 43 does not characterize uniquely f z Condition 2 is more difficult to analyze because it is not convex In order to solve this problem numerically we approx imate condition 2 with a convex constraint Condition 2 defines the value of the wavelet transform at the points x3 iez Instead of imposing that the local maxima of W h x are located at these points we impose that W gt A x is as small as possible on average This generally creates local modulus maxima at the positions x n z2 The number of modulus maxima of Ws f x depends on how much this function oscillates To have as few modulus maxima as possible outside the abscissa x j n ez2 we also minimize the energy of the derivative of W h x Since these conditions must be imposed at all scales 27 we minimize globally 2 Wal Weih e e2 00 dWh gt wane P an j The weight 2 expresses that the relative smoothness of W f x increases with the scale 27 Let t x be the derivative of x If there exist two constants A gt gt 0 and B such that for all w R 00 00 Ar lt So Wo S gt Ww lt Bo 48 j then for any h x L R Ag h lt IAI lt BallAll j x 49 Hence is a norm over L
6. 1 3 In order to have a wavelet antisymmetrical with respect to 0 and z symmetrical with respect to 0 the shifting constant w of 92 is equal to 1 2 Equations 101 prove that w x is a quadratic spline with compact support whereas a is a cubic spline whose integral is equal to 1 The 27 periodic function H w G w and K w can be viewed as the transfer function of discrete filters with finite impulse response The corresponding impulse responses are given in Table I These filters are used in fast wavelet transform computations APPENDIX B FAST WAVELET ALGORITHMS FOR 1 D SIGNALS This appendix describes an algorithm for computing a discrete wavelet transform and the inverse algorithm that reconstructs the original signal from its wavelet transform We suppose that the wavelet x is characterized by the three discrete filters H G and K described in Appendix A We denote H G and K the discrete filters obtained by putting 2P 1 zeros between each of the coefficients of the filters H G and K The transfer function of these filters is respectively H 2 w G 2 w and K 2 w We also denote by Hp the filter whose transfer function is the complex conjugates of H 2 w H 2 w We denote by A x B the convolution of two discrete signals A and B MALLAT AND ZHONG CHARACTERIZATION OF SIGNALS The following algorithm computes the discrete wavelet transform of the discrete signal S f At each scale 2 it decomposes S4
7. Equation 87 proves that V C K The sequences of functions that satisfy condition 1 are the elements of K that belong to A VNTI To reconstruct the element of IT N V that minimizes the norm we alternate projections on I and V that are orthogonal with respect to the norm As in one dimension one can prove that the orthogonal projection on V is the operator Py WoW that was defined by 64 The orthogonal projections P p on T are defined in Appendix E For a discrete image of N pixels the implementations of both Py and Pp require O N logs N operations Let P PyoPp be the alternate projection on both sets Since I is an affine space and V a vector space for any initial sequence X 9 2 4 95 Y ez lim P X converges strongly to nm 00 725 70 96 60 25 4 50 2 98 40 92 30 21 20 10 0 20 40 60 80 100 a 40 35 30 25 20 15 10 O 20 40 60 80 100 b Fig 10 a SNR when reconstructing the wavelet transform components W LF and We i f from the modulus maxima of the Lena image shown in Fig 9 The abscissa gives the number of iterations on the operator P Each curve is labeled by the scale 2 for 1 lt j lt 6 b SNR of the reconstructed Lena images computed with respect to the original image as a function of the number of iterations on the operator P the orthogonal projection of X onto A INV Hence if we begin the iteration from the zero element of K the algorithm converges strongly
8. it is important to introduce no distortion around the eyes because these are highly visible for a human observer In the following we do not introduce such context information in the selection To code efficiently the edge information we need to take ad vantage of the similarities between edges obtained at different scales As it can be observed in Fig 9 the edges of the main image structures have similar positions at the three finer scales 91 92 and 2 These three finer scales also carry more than 90 of the image frequency bandwidth and thus cover most of the image information We build our edge encoding from these scales only The coarse scale information corresponding to the wavelet transform at scales 2 gt 2 is kept as a low frequency image Se f defined in Section VI B The edge selection is first performed at the scale 2 because at the finer scale 21 edges are too greatly contaminated by high frequency noises The boundaries of the important coherent structures of ten generate long edge curves We thus remove any edge curve whose length is smaller than a given length threshold Among the remaining curves we select the ones that correspond to the sharpest image variations This is done by removing the edge curves along which the average value of the wavelet transform modulus is smaller than a given amplitude threshold MALLAT AND ZHONG CHARACTERIZATION OF SIGNALS Once the selection is done we must efficiently co
9. scales and detect sharp variation points from their first or second order derivative The extrema of the first derivative correspond to the zero crossings of the second derivative and to the inflection points of the smoothed signal This section explains how these multiscale edge detection algorithms are related to the wavelet transform We call a smoothing function any function x whose integral is equal to 1 and that converges to 0 at infinity For example one can choose equal to a Gaussian We suppose that x is twice differentiable and define respectively Y x and 7 a as the first and second order derivative of 0 z dO x B d 6 a w x aa and Y a J72 1 711 By definition the functions x and y x can be consid ered to be wavelets because their integral is equal to 0 00 00 we a dz 0 and y x dx 0 0 oO In this paper we denote which is the dilation by a scaling factor s of any function x A wavelet transform is computed by convolving the signal with a dilated wavelet The wavelet transform of f z at the scale s and position x computed with respect to the wavelet w x is defined by Wee i eo 2 The wavelet transform of f x with respect to w x is We f x f p z 3 We derive that We f x fx so d d Na Ff 8 e and 4 T 2 2 Tyla EUa The wavelet transforms W2 f x and W f x are respec tively the first and second deri
10. therefore the parameters K o and g do not change much We thus estimate these values for portions of chains by looking at the evolution of the modulus values across scales Let us suppose that we have a portion of the maxima chain that propagates between the scales 2 and 2 We also suppose 724 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE VOL 14 NO 7 JULY 1992 that in a given neighborhood at each scale 27 the value of Mo f x y is bounded by its values along this maxima chain This means that the maxima chain corresponds to the sharpest image variation in the neighborhood Since Mo f x y is bounded by the maxima values we estimate the parameters a and o that satisfy 80 from the evolution across scales of these modulus maxima values In theory this should be done by using the absolute maximum of M gt f x y along the maxima chain at each scale 27 It is often better to regularize these computations by averaging the maxima modulus value along the corresponding portion of the chain This is justified since we suppose that the singularity type does not vary greatly along this portion of chain Let a be the average value of Mo f x y As in one dimension we estimate the smoothing factor g and the Lipschitz regularity ag by computing the values that minimize I gt iog ateni j l ag l 2 logs a 2 81 This algorithm associates with each portion of maxima chain three constants K
11. we prove this implication by applying the Parseval theorem on each L R norm component of the norm defined in 85 We saw that the set of functions A z y whose wavelet transform satisfies condition 1 is the closed affine space f O The minimization of the norm over this closed convex has a unique solution whose computation might however not be stable As in one dimension we can prove that the computation of this minimum is stable if and only if the family of functions T apa Ge T yi _ y T a t yi g Wb Geezer 88 is a frame of the space U that they generate The factor 2 normalizes the L R norm of the functions in this family The frame condition expresses the equivalence of the L R norm of any function in U and the sum square of the inner products of this function with each function of the family 88 To compute the solution of our minimization problem we use an alternate projection algorithm just as in one di mension Let K the space of all sequences of the function 9 2 y 93 2 y such that 9 y 92 2 9 en lt 00 where the norm is defined by expression 85 We define the set I of all sequences of functions 95 2 9 95 7Y jez K such that for any index j and all maxima position x7 y 95 al y Wo f z3 yt and g x y2 W3 f 21 y2 89 The set I is an affine space that is closed in K Let V be the space of dyadic wavelet transforms of all functions in L R
12. Each curve is labeled by the scale 2 for 1 lt lt 6 b SNR of the reconstructed signal computed with respect to the original signal as a function of the number of iterations on the operator P by increasing the number of iterations is negligible Since each iteration requires O N log N operations these recon structions do not require extensive computations and can be done in real time The reconstructed functions are not equal to the original signal but are numerically close They have no spurious oscillations and the same types of sharp variations Qualitatively the reconstructed signals are thus very similar to the original one and the errors are hardly noticeable by comparing the graphs as shown by Fig 5 We have no upper bound on the error due to the distance between the signal to which we converge and the original signal This is an open mathematical problem but the numerical precision of this reconstruction algorithm is sufficient for many signal processing applications VI WAVELET TRANSFORM OF IMAGES We explained in Section H that in two dimensions a mul tiscale edge detection can be reformalized through a wavelet transform defined with respect to two wavelets y x y and y x y The second part of this article extends our 1 D results for image processing applications A General Properties We denote that y x y sv 3 5 and that p er y w5 4 The wavelet transform of a function f x y L R
13. R which is equivalent to the classical L7 R norm We prove that 48 implies 49 by observing that dWoas h a a zM Sf rina As for the norm equivalence equation 15 we then prove the implication by applying the Parseval theorem to each L R norm component of the norm defined by 47 Equation 48 is valid for any dyadic wavelet y x that satisfies 12 and is smooth enough For example there exist two such constant A and Bo for the wavelet shown in Fig 1 a By replacing condi tion 2 by the minimization of A we define a problem that has a unique solution Indeed condition 1 imposes that h a must belong to the closed affine space f O and the minimization of a norm over such a closed convex has a unique solution Although there exists a unique element of f O whose norm is minimum the computation of this function might not be stable If the two constants Aj and B of 48 are equal the norm is proportional to the classical L R norm The solution of the minimization problem is therefore the orthogonal projection of f x over U The frame theory proves 6 that one can make a stable computation of the orthogonal projection of f x onto U from the inner products f x pa x3 6 nyez if and only if the family of functions Va a 2 o is a frame of U The JNJEL 50 factor V2 normalizes the L R norm of each function in the family By definition such a family is a frame of U 6 if a
14. a review by Hummel and Moniot 11 for more details If the smoothing function x is a Gaussian the properties of the wavelet transform zero crossings are more easily un derstood because W f x can be interpreted as the solution of a heat diffusion process at time t s 12 With this approach Hummel and Moniot 11 as well as Yuille and Poggio 27 have proved some completeness properties under restrictive conditions like supposing that f x is a polynomial MALLAT AND ZHONG CHARACTERIZATION OF SIGNALS 29 In general there are known counterexamples that prove that the positions of the zero crossings of W2 f x do not characterize uniquely the function f a For example the wavelet transforms of sin x and sin x 2 sin 22 have the Same zero crossings at all scales s gt 0 Meyer 21 found a large family of such counterexamples To obtain a complete and stable zero crossing representa tion Mallat 16 conjectured that it is sufficient to record the zero crossing positions of W2 f x at all dyadic scales 2 f as well as the integral values intl C We f u du 39 between any pair of consecutive zero crossings Zn Zn 1 This conjecture was motivated by a reconstruction algorithm that is able to recover a close approximation of the original signal from these zero crossings and integral values 18 We proved in Section II that the zero crossings of W3 f x occur at the extrema points of the wavelet transform W5 f
15. any scale s gt 0 Hence ao is the upper bound of the set of a such that there exists K that satisfy IWsh x lt Ks 716 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE VOL 14 NO 7 JULY 1992 for any scale s gt 0 and any x in the corresponding neighbor hood of xo By inserting this inequality in 35 we obtain Wa f a lt K2 s07 with so V233 072 37 This equation is satisfied at all points if and only if it is satisfied at the locations of all the local maxima of W f x If the signal is multiplied by a constant then K is also multiplied by but c and ap are not affected On the contrary if the signal is smoothed by a Gaussian of variance o and integral 1 then AK and ap are not affected but c becomes a a4 This shows clearly that the parameters a o and K describe different properties of the sharp variation that occurs at xo Fig 3 gives the examples of a step edge and a Dirac smoothed by Gaussians of different variances The decay of the maxima are clearly affected by the different Lipschitz exponents as well as the variance of the Gaussian smoothing Let us explain how to compute numerically the Lipschitz regularity and the smoothing scale o from the evolution of the wavelet transform modulus maxima across scales If we detect the modulus maxima at all scales s instead of just dyadic scales 2 their position would define a smooth curve in the scale space plane s 2 These c
16. ao and o that describe the intensity profile of the image sharp variation along the chain Such a characterization of edge types is important for pattern recognition For example we can discriminate occlusions from shadows by looking at whether the image intensity is discontinuous or is smoothly varying For the circle image of Fig 8 the wavelet transform modulus along the boundary remains constant across scales which means that ag QO and 0 Indeed the image intensity is discontinuous along the border and the constant K gives the amplitude of the discontinuity In general we believe that an edge detection should not be viewed as a binary process that labels the image pixels as edge points or nonedge points but as a procedure that characterizes precisely the different types of image sharp variations VIII RECONSTRUCTION OF IMAGES FROM MULTISCALE EDGES A Reconstruction Algorithm The algorithm that reconstructs images from the local max ima of their wavelet transform modulus is an extension of the 1 D algorithm described in Section V B Let f x y L R and W4 f x y W2 f a Y sez be its dyadic wavelet trans form For each scale 27 we detect the local maxima of Mo f x y along the direction given by the angle image A f x y We record the positions of the modulus maxima elvd op as well as Ma f x y Az Fi vi ep In two dimensions the number of modulus maxima is no longer countable From M f x yi
17. f into Sz 1 f and Wz 1 f 7 20 while j lt J Woy aod i Sof Gj Spd z S f H p zt end of while The proof of this algorithm is based on the properties of the wavelet y x described in Appendix A At each scale 2 we divide the values of the samples of 8 f Gj by j to obtain accurate measures of Lipschitz exponents from the wavelet maxima Due to discretization the wavelet modulus maxima of a step edge do not have the same amplitude at all scales as they should in a continuous model The constants A compensate for this discrete effect The values of that correspond to the filters of Table I are given in Table II The border problems are treated by making symmetry and a periodization of S1 f n lt lt y as explained in Section III B The convolutions must take into account this periodization The complexity of this discrete wavelet transform algorithm is O N log N and the complexity constant is proportional to the number of nonzero coefficients in the impulse response of the filters H and G The inverse wavelet transform algorithm reconstructs S f from the discrete dyadic wavelet transform At each scale 25 it reconstructs S4 _ f from S4 f and W4 f The complexity of this reconstruction algorithm is also O N log N Td while j gt 0 Che j W f Kj 1 Saf Hiz i io age end of while APPENDIX C A PARTICULAR CLASS OF 2 D DYADIC WAVELETS In this appendix we characterize th
18. how to characterize different types of sharp variation points from the evolution across scales of W f x at the modulus maxima locations The Canny edge detector is easily extended in two dimen sions We denote by We f a f 8 ree sials y EG 712 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE VOL 14 NO 7 JULY 1992 the dilation by s of any 2 D function x y We use the term 2 D smoothing function to describe any function 6 z y whose integral over z and y is equal to 1 and converges to 0 at infinity The image f x y is smoothed at different scales s by a convolution with 6 x y We then compute the gradient vector V f 6 x y The direction of the gradient vector at a point Zo yo indicates the direction in the image plane x y along which the directional derivative of f z y has the largest absolute value Edges are defined as points o yo where the modulus of the gradient vector is maximum in the direction towards which the gradient vector points in the image plane Edge points are inflection points of the surface f 0 x y Let us relate this edge detection to a 2 D wavelet transform We define two wavelet functions 7 x y and w x y such that 1 _ 06 x y 2 _ 00 z y v e B and P e y Let pi z y aw 2 and 2 z y 4y 2 4 Let f z y L R The wavelet transform of f x y at the scale s has two components defined by W3 f z y fxb3 a y and W f
19. maxima of nonisolated singularities is studied in more detail in 17 V SIGNAL RECONSTRUCTION FROM MULTISCALE EDGES Section IV shows that one can get a precise description of the signal sharp variation points from the evolution of the wavelet transform modulus maxima across scales An important question is to understand whether the whole signal information is embedded into these modulus maxima Is it possible to have a stable signal reconstruction only from the modulus maxima information at the dyadic scales 27 The next section reviews briefly some results on the reconstruction of signals from zero crossings and multiscale edges Section V B describes an algorithm that reconstructs a close approximation of the original signal from the wavelet transform modulus maxima Numerical results are presented in Section V C A Previous Results The reconstruction of signals from multiscale edges has mainly been studied in the zero crossing framework We saw in Section II that if the wavelet is given by w x aly multiscale edges are detected from the zero crossings of the wavelet transform W f a The most basic result concerning the reconstruction of signals from the zero crossing is the Logan theorem 14 However as it is explained in 18 the hypotheses of the Logan theorem are not appropriate to study the reconstruction of signals from multiscale edges The Logan theorem has been generalized by several authors 4 24 28 Refer to
20. of iterations If the original signal has NV samples we record the positions and values of the modulus maxima at all scales 2 for 1 lt j lt log NV 1 We also keep the average value of the original discrete signal which characterizes S f for J log N 1 as explained in Section II B Equation 53 proves that we can compute the projection Py by implement ing W followed by W With the fast algorithms described in Appendix B this requires a total of O N log N operations Appendix E proves that the implementation of Pp also requires O N logo N operations The projection operator that suppresses the wavelet transform oscillations is computed with the same complexity Hence each iteration on P involves O N logo N operations The rms signal to noise ratio SNR of the reconstruction is measured in decibels At the scale 27 for 1 lt j lt 6 Fig 6 a gives the value of the SNR for the reconstruction of W f after n iterations on the operator P with 1 lt n lt 100 At all scales the error decreases quickly during the first 20 iterations and then decays much more slowly For a fixed number of iterations on P the SNR increases when the scale increases This proves that the remaining error is rather concentrated at fine scales like in the counterexamples of Meyer 21 After n iterations we can reconstruct a signal by applying the inverse wavelet transform operator on the reconstructed wavelet transform Fig 6 b sho
21. replace the operators W and W by their expression given in 9 and 16 we obtain 00 VIEL g K g z with 18 l oo K Xai Pai z 19 These equations are known as reproducing kernel equations The energy of the kernel K x measures the redundancy of the wavelet transform at the scales 2 and 2 MALLAT AND ZHONG CHARACTERIZATION OF SIGNALS 0 8 1 4 0 0 0 8 0 0 2 1 2 42 1 0 1 2 a b Fig 1 a This wavelet is a quadratic spline of compact support that is continuously differentiable It is defined in Appendix A and it is the derivative of the cubic spline function 6 x shown in b Fig 1 a is a quadratic spline wavelet of compact support which is further defined in Appendix A It is the derivative of the smoothing function x shown in Fig 1 b Fig 2 a is the plot of a discrete signal of 256 samples Fig 2 b shows its discrete dyadic wavelet transform computed on nine scales At each scale 2 we compute a uniform sampling of the wavelet transform that we denote we f The next section explains how to discretize the continuous wavelet model and solve border problems Fast algorithms to compute the wavelet and the inverse wavelet transform are described in Appendix B The reader not interested in numerical issues might want to skip Section III B Since our wavelet is the derivative of a smoothing function 4 proves that we f is proportional to the derivative of the original s
22. to the element of A whose norm is minimum B Numerical Reconstruction of Images from Multiscale Edges We study the error of the reconstruction algorithm as a function of the number of iterations on the operator P At each scale 27 the SNR integrates the error on the horizontal and the vertical components of the wavelet transform Fig 10 a gives the evolution of the SNR when reconstructing the wavelet transform of the Lena image from the modulus maxima shown in the Fig 9 After n iterations we reconstruct an image by applying the inverse wavelet transform operator on the reconstructed wavelet transform Fig 10 b is the SNR of the reconstructed images computed with respect to the original one as a function of the number of iterations on the operator P The graphs of Fig 10 are very similar to the graphs of Fig 6 that show the reconstruction SNR for a 1 D signal The increase is fast during the first 20 iterations and then slows down After a given number of iterations Fig 10 a shows that the error is mostly concentrated at fine scales This error has two components The first one is the distance to the wavelet transform to which we converge and the other one is the distance between the point to which we converge and the wavelet transform of the original image As in one dimension the convergence is exponential but the convergence rate is very slow After 20 iterations the distance between the reconstructed image and the imag
23. y and bley given in Appendix C from the sample values 5 fn Men mez we can compute a uniform sampling of the wavelet transform of f z y at scales larger than 1 For any scale 2 we denote i yf Fa p Pi H pa j H a fi rL pi FPL c ni S s i x an r r i da a 4g ba i ETE g P oA mr P Wee f We fin a mm e t AR ge T gpd AE vale nome d ered Fig g gs Ps j i i Sk OS F om S85 fla d u A wy seas 22 PO tae PUP wW mop N fiin mice The sampling shift w depends on the choice of wavelets Given an image D Appendix D describes a fast algorithm MALLAT AND ZHONG CHARACTERIZATION OF SIGNALS that computes the discrete dyadic wavelet transform Shes Wo F ey Co o 72 Images are finite 2 D discrete signals D dnm n m EN of N by N pixels We solve border problems such as those in a 2 D cosine transform We suppose that the image is symmetrical with respect to each of its border and has a period of 2N by 2N pixels For J log N 1 one can prove that So z iS constant and equal to the average of the original image D Figs 8 and 9 give two examples of wavelet transforms computed over log V 1 scales The numerical complexity of the fast discrete wavelet transform is O N log N The reconstruction of the original image from its discrete wavelet transform is also performed with O N log N operations The discrete modulus images M4 f and angle images Ad f are computed with 67 and 68 The m
24. z y Similarly to 4 one can easily prove that tad 3 f s z Ta miren 3 Ep o e y ZVO A Hence edge points can be located from the two components Wi f x y and W2 f z y of the wavelet transform fay a y 7 II DYADIC WAVELETTRANSFORM IN ONE DIMENSION A General Properties For most purposes the wavelet model is not required to keep a continuous scale parameter s To allow fast numerical implementations we impose that the scale varies only along the dyadic sequence 27a jez We review the main properties of a dyadic wavelet transform and explain under what con dition it is complete and stable For thorough presentations of the wavelet transform refer to the mathematical books of Meyer 20 and Daubechies 5 or to signal processing oriented reviews 15 22 The wavelet model has first been formalized by Grossmann and Morlet 10 A wavelet is a function z whose average is zero We denote by wo x the dilation of y x by a factor 2 ya 2 ZWS The wavelet transform of f x at the scale 2 and at the position x is defined by the convolution product Wo f x f Yz a 9 We refer to the dyadic wavelet transform as the sequence of functions Wf Wz f 2 jez and W is the dyadic wavelet transform operator 10 Let us study the completeness and stability of a dyadic wavelet transform The Fourier transform of W3 f x is Wo f w f oj w 11 By imposing that there exists
25. 710 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE VOL 14 NO 7 JULY 1992 Characterization of Signals from Multiscale Edges Stephane Mallat and Sifen Zhong Abstract A multiscale Canny edge detection is equivalent to finding the local maxima of a wavelet transform We study the properties of multiscale edges through the wavelet theory For pattern recognition one often needs to discriminate different types of edges We show that the evolution of wavelet local maxima across scales characterize the local shape of irregular structures Numerical descriptors of edge types are derived The completeness of a multiscale edge representation is also studied We describe an algorithm that reconstructs a close approximation of 1 D and 2 D signals from their multiscale edges For images the reconstruction errors are below our visual sensitivity As an application we implement a compact image coding algorithm that selects important edges and compresses the image data by factors over 30 Index Terms Edge detection feature extraction level cross ings multiscale wavelets I INTRODUCTION OINTS OF SHARP variations are often among the most important features for analyzing the properties of transient signals or images In images they are generally located at the boundaries of important image structures In order to detect the contours of small structures as well as the boundaries of larger objects several researchers in comp
26. AND MACHINE INTELLIGENCE VOL 14 NO 7 JULY 1992 that Va zi x that h x g x I is a frame of U 53 implies nj EZ 1 eren Ae S gt 2 A x glr va ad e l n j EZ 128 Since Wo A 2 ez E I we have W h x Wo f xi and thus h x 9 2 boi 23 2 Wo f 24 Wa g a a 129 From the norm equivalence of 49 we can also derive that 2 Wahle jez Wala jez lt Ballh x ge 130 Since the projector P 4 is orthogonal and Wa hA x ez E A 2 Wagt ejz PAWagl jezl 2 lt Wot jez E Wa h x jez 131 Equations 128 131 imply that 2 Wai jez PalWa9 je2 Bo Pre lt 2 E Poe 132 A3 ene n j EZ This proves 120 for Cp 42 From 123 and 132 we then derive that DA 2B This inequality gives a lower bound for the angle between the affine space I and the space V Let P PyoPp be the alternate projection on both spaces A classical result on alternate projections 26 enables us to derive that for any element X K there exists a constant R such that DA Pix SPO x Shi A E IV X Ppx IX P X 133 Pie 134 This proves that the algorithm converges exponentially with a D 1 2 convergence rate larger than 1 55 gt a REFERENCES 1 Z Berman The uniqueness question of discrete wavelet maxima representation Tech Rep TR 91 48 Syst Res Cent Univ of
27. INTELLIGENCE VOL 14 NO 7 JULY 1992 described a compact image coding procedure that selects the important visual information before coding The compression rates are between 30 and 40 in the examples that are shown but most of the light image textures are not coded A double layer coding based on multiscale edges and textures has recently been developed by Froment and Mallat 8 APPENDIX A A PARTICULAR CLASS OF 1 D WAVELETS This appendix defines a class of wavelets that can be used for a fast implementation of discrete algorithms We first define the smoothing function z introduced in Section III B and then we build the wavelet y x and the reconstructing wavelet x x which are associated with 2 We impose that the Fourier transform of the smoothing function x defined by 20 can be written as an infinite product 00 dw e H 2 p 1 90 where H w is a 27 periodic differentiable function such that H w H w 7 lt 1 and A 0 1 91 One can prove 15 that the conditions 91 are sufficient so that 90 defines a smoothing function x which is in L R The parameter w is the sampling shift that was introduced in Section III B It is adjusted in order so that 2 is symmetrical with respect to 0 Equation 90 implies that 2w e H w p w We define a wavelet y x whose Fourier transform 7 w is given by 92 P 2w e G w d w where G w is a 27 periodic fun
28. Maryland College Park Apr 1991 2 J Canny A computational approach to edge detection JEEE Trans Patt Anal Machine Intell vol PAMI 8 pp 679 698 1986 3 S Carlsson Sketch based coding of grey level images Signal Pro cessing vol 15 pp 57 83 July 1988 4 S Curtis S Shitz and V Oppenheim Reconstruction of nonperiodic two dimensional signals from zero crossings IEEE Trans Acoustic Speech Signal Processing vol 35 pp 890 893 1987 5 I Daubechies Ten lectures on wavelets CBMS NSF Series Appl Math SIAM 1991 6 R Duffin and A Schaeffer A class of nonharmonic Fourier series Trans Amer Math Soc vol 72 pp 341 366 1952 7 G Folland Introduction to Partial Differential Equations NJ Princeton University Press 1976 8 J Froment and S Mallat Second generation compact image coding with wavelets in Wavelets A Tutorial in Theory and Applications C Chui Ed New York Academic Jan 1992 pp 655 678 Princeton 9 K Grochenig Sharp results on random sampling of band limited functions in Proc NATO ASI Workshop Stochastic Processes New York Kluwer Aug 1991 A Grossmann and J Morlet Decomposition of Hardy functions into square integrable wavelets of constant shape SIAM J Math vol 15 pp 723 736 1984 R Hummel and R Moniot Reconstruction from zero crossings in scale space JEEE Trans Acoustic Sp
29. after projections on I and V Hence the operator Py acts as the identity operator and P can be rewritten P PyoP p The analysis of Section V B proves that we are then guaranteed to converge strongly to an element in Y AT AV In two dimensions the operator P p transforms a sequence g5 2 9 97 29 seg K into the closest sequence hilz y hlz y ez E T Let ej x y 5 2 9 ez MALLAT AND ZHONG CHARACTERIZATION OF SIGNALS be such that for any j Z amp z y gj z y hj x y and x y g x y A x y The sequence W x y h E seg is chosen so that 00 2 CRIE He de 114 is minimum The constraints on e x y and 7 x y are inde pendent The minimization of 114 is obtained by minimizing each component 142 lhe IP yji p 115 and de2 2 2 4 923 2 l l 22115 zal for all integers 7 Z Let us concentrate on the minimization of 115 Let xo y and 21 y be two consecutive modulus maxima position at a fixed y The function z y must satisfy l Zo y e 1 y 116 Wo f xo y g 0 y Wi f 21 y gj 1 9 The minimization of 115 subject to these constraints is obtained by minimizing J i To For y fixed we obtain a 1 D minimization problem which is identical to the minimization of expression 111 The solution is a sum of two exponentials as in 113 This analysis shows that the solution of the 2 D minimization probl
30. alues This chaining procedure defines an image representation that is a set of maxima chains Image edges might correspond to very different types of sharp varia tions As in one dimension we discriminate different types of singularities by measuring their local Lipschitz regularity Definition 2 Let 0 lt a lt 1 A function f x y is uniformly Lipschitz over an open set 2 of R if and only if there exists a constant K such that for all xo yo and 1 Y1 in Q f 0 yo f z1 y1 lt K xo 21 yo 1 73 The Lipschitz regularity of f x y over Q is the superior bound of all such that f x y is uniformly Lipschitz a 723 In two dimensions the Lipschitz regularity is character ized by the decay across scales of both W f x y and W f x y The decay of these two components is bounded by the decay of Mp f x y Let us suppose that the two wavelets y x y and a y are continuously differentiable and that their decay at infinity is Olapar Theorem 2 Let 0 lt a lt 1 A function f x y is uniformly Lipschitz over an open set of R if and only if there exists a constant K such that for all points x y of this open set Mo f x y lt K 27 74 This theorem is the 2 D extension of Theorem 1 and its proof is essentially the same 20 The logarithm of 74 yields loga M f x y lt loga K aj Uniform Lipschitz exponents can thus be measured from the evolution a
31. and A f x y1 we can compute W4 f 23 43 W2 f x2 y2 and vice versa The inverse problem consists of finding the set of functions h x y that satisfy the following two constraints 1 At each scale 2 and for each modulus maxima location zi y we have Wl h x yi Wy 21 y2 and W3h xi y2 WZ f 23 y2 2 At each scale 27 the modulus maxima obtained from Wy h x y and W3 h x y are located at the abscissa i yi vER Let us analyze property 1 At any point x9 yo the wavelet transform can be rewritten as inner products W3 h xo yo f z y pl xo T Yo y W2 h z0 yo f z y pz xo T Yo y Let U be the closure of the set of functions that are linear combinations of any function of the family 82 25 a5 x aj z yi 7 y pz E T yf g WS GeyezxR The factor 2 normalizes the L R norm of each function One can prove that the set of functions A x y whose wavelet transform satisfies the condition 1 are the functions whose or thogonal projection on U is equal to the orthogonal projection of f x y on U Let O be the orthogonal complement of U in L R This set is therefore the affine space f O of functions that can be written h x y f x y g x y with g z y O 84 We replace condition 2 with a convex constraint that has a similar effect in order to solve the problem numerically We do not impose that the points 23 y3 Gojezxr ate the on
32. as f xz 2 y Wo f X2 2 Y J 62 A 2 D dyadic wavelet transform is more than complete it is redundant Any sequence of functions g x y g 2 WD se 5 is not necessarily the dyadic wavelet transform of some functions in L R We denote by W the operator defined by Ww 9 2 y 9 z y ez 2 lt Ba 721 o E g xz 2 y 97 x55 ay 63 J The sequence 9 2 4 95 2 9 cz is a dyadic wavelet transform if and only if w wo 93 x y 95 2 sez 95 2 9 95 Y jez 64 In Section I we explained that multiscale sharp variation points can be obtained from a dyadic wavelet transform if OO x y OO x y Ox l p z y and 4 x y ae 65 Equation 8 proves that the wavelet transform can be rewritten rie Eire UV S 09 a y 66 The two components of the wavelet transform are proportional to the two components of the gradient vector V f 65 2 y This appears clearly in Fig 8 which shows the 2 D wavelet transform of the image of a circle At each scale 2 the modulus of the gradient vector is proportional to Mo f x y The angle of the gradient vector with the horizontal direction is given by Ag f a y Like in the Canny algorithm 2 the sharp variation points of f x f z y are the points x y where the modulus Mo f x y has a local maxima in the direction of the gradient given by A f x y We record the position o
33. ations in signal processing It allows us to process the image information with edge based algorithms We describe a compact image coding algorithm that keeps only the important edges The image that is recovered from these main features has lost some small details but is visually of good quality Examples with compression ratio over 30 are shown Another application to the removal of noises from signals is described in 17 The article is organized as follows Section II relates mul tiscale edge detection to the wavelet transform It shows that a Canny edge detector is equivalent to finding the local maxima of a wavelet transform modulus Until Section VI we concentrate on 1 D signals Section III A reviews the wavelet transform properties that are important for under standing multiscale edges The wavelet transform is first defined over functions of continuous variables and Section H B explains how to discretize this model The numerical implementation of fast wavelet transform algorithms is given in Appendix B Section IV explains how to characterize different types of sharp signal variations from the evolution across scales of the wavelet transform maxima Section V studies the reconstruction of signals from multiscale edges We review some previous results and explain how to formalize the reconstruction problem within the wavelet framework The reconstruction algorithm is described in Section V B and numerical results are presented i
34. avelet maxima across scales Definition 1 Let 0 lt a lt 1 A function f z is uniformly Lipschitz over an interval a b if and only if there exists a constant K such that for any zo 1 Ja b f zo f a1 lt K zo x1 We refer to the Lipschitz uniform regularity of f a as the up per bound ao of all such that f x is uniformly Lipschitz a If f x is differentiable at xo then it is Lipschitz a 1 If the uniform Lipschitz regularity apo is larger the singularity at r will be more regular If f x is discontinuous but bounded in the neighborhood of zo its uniform Lipschitz regularity in the neighborhood of xo is 0 Theorem 1 proves that the Lipschitz exponent of a function can be measured from the evolution across scales of the absolute value of the wavelet transform We suppose that the wavelet x is continuously differentiable and has a decay at infinity that is O Theorem 1 Let 0 lt a lt 1 A function f z is uniformly Lipschitz over Ja b if and only if there exists a constant K gt 0 such that for all z a b the wavelet transform satisfies 27 Wz f x lt K 27 28 The proof of this theorem can be found in 20 From 28 we derive that logs Wo f x lt logo K a3 If the uniform Lipschitz regularity is positive 28 implies that the amplitude of the wavelet transform modulus maxima should decrease when the scale decreases On the contrary the
35. cross scales of log Mo f x y This result enables us to discriminate between different types of singu larities When the signal variations are smooth we can measure how smooth they are with the same approach as in one dimension Locally we model the smooth variation of f x y at xo yo as the convolution of a function h x y that has a singularity at o Yo with a 2 D rotationally symmetric Gaussian of variance g 75 a y Fo 76 We suppose that the uniform Lipschitz regularity of h x y in a neighborhood of zo yo is ag If the two wavelets Y x y and x y are the partial derivatives of a smoothing function 6 x y which closely approximates a rotationally symmetric Gaussian then we can estimate the variance o The wavelet transform modulus of f x y is defined at any scale s by Ms f x y yIW f x y W2 f z y With the same derivations as for 35 we prove that 1 Jr y h go T y with Gly Zo exp 77 J Mo f x y amp Mehle with so 227 07 78 0 Equation 74 of Theorem 2 is valid not only at dyadic scales 2 but at all scales s gt 0 For lt ag h x y is uniformly Lipschitz in a neighborhood of xo yo Hence there exists K gt 0 such that for any points x y in this neighborhood M h z y lt Ksg 19 We thus derive from 78 that Mo f x y lt K2 s871 with so V2 07 80 Along a maxima chain the singularity type varies smoothly
36. ction Equation 22 for J 1 proves that w w x w and w must satisfy 93 th 2w X 2w llw o 2w 94 Let us impose that w can be written X 2w e K w d w 95 where K w is a 27 periodic function If we insert 92 and 95 into 94 we obtain H w G w K w 1 96 One can prove that 96 is sufficient to define K w such that X w is the Fourier transform of a reconstructing wavelet that satisfies 13 We want a wavelet y x equal to the first order derivative of a smoothing function x This implies that 7 u must have a zero of order 1 at w 0 Since 4 0 1 93 yields that G w must have a zero of order 1 at w 0 We choose H w in order to obtain a wavelet x which is antisymmetrical as regular as possible and has a small compact support A TABLE I FINITE IMPULSE RESPONSE OF THE FILTERS H G N AND L THAT CORRESPOND TO THE QUADRATIC SPLINE WAVELET OF Fic 1 a 0 0078 125 0 054685 0 171875 0 0078125 0 046875 0 1171875 0 65625 0 1171875 0 046875 0 0078125 0 171875 0 054685 0 0078125 family of 27 periodic functions H w G w and K w that satisfy these constraints is given by H w e cos w 2 97 G w 4ie sin w 2 98 H 2 K w a 99 From 90 and 93 we derive 2n 1 bw ee 2 100 2n 2 p w Se 101 The Fourier transform 6 w of the primitive is therefore i 2n 2 6 w ae 102 In the example of Fig 1 we chose 2n
37. de the remaining information This requires coding the position modulus value and angle value of each modulus maximum along the maxima curves at the scales 21 2 and 2 plus the low frequency image SS f The geometry of the edge curves is coded only at the scale 2 because we neglect the differences between the maxima positions at the scales 2 2 and 2 Maxima chains are coded by recording the position of the first point of cach chain and then coding the mcrement between the position of one edge point to the next one along the chain Carlsson 3 showed that this requires on average 10 b for the first point and 1 3 b per point along the chain with an entropy coding At each scale the direction of the gradient image intensity at the edge locations is approximately orthogonal to the tangent of the edge curves We thus do not code the angle values but approximate each of them by the orthogonal direction of the edge tangent at the corresponding location The values of the modulus along the edge curves at the scales 2 27 and 2 are recorded with a simple predictive coding using a coarse quantization of the prediction values In the frequency domain the image S f has an energy that is mostly concentrated in a domain that is 2 times smaller than the frequency support of the original image We thus subsample this image along its rows and columns by a factor 2 and its grey level values are coded on 6 b We give in Fig 12 two e
38. e 2 D wavelets used for numerical computations In order to compute the wavelet transform with a minimum amount of operations we choose two wavelets w x y and w z y that can be written as separable products of functions of the x and y variables Let y x be a wavelet that belongs to the class described in Appendix A and whose Fourier transform is defined by 2w e G w o w with do w 00 eww I HO Pw 103 p i We define Y z y p x 2bo 2y and Y z y 2o 2x h y 104 Since rz ole these two wavelets can be rewritten pir y Zw 007 x y f 2 j and w x y with 105 729 TABLE II NORMALIZATION COEFFICIENTS A FOR THE QUADRATIC WAVELET OF Fic 1 a For 7 gt 5 A 1 6 x y O x 2ho 2y and 07 x y 260 2x 0 y The fast discrete wavelet algorithm does not allow us to have 61 a y 07 x y However the functions 61 z y and 0 x y are numerically close enough so that they can be considered to be equal to a single function 6 x y in a first approximation Equations 103 and 104 imply that the Fourier transform of the two wavelets w1 a y and Y x y are given by Y 2wa 2wy 67 Gwe bo we Go wy and y 2wr 2wy bo wa e vs G wy do wy One can define the smoothing function z y of Section VI B as 106 p z y bo x dbo y Let us now define wz wy and wz wy such that R 2wa 2wy eM K We L wy we wy
39. e Wa f x is obtained through a convolution with 72 2 it is also differentiable in the sense of Sobolev and it has at most a countable number of modulus maxima Let x nez be the abscissa where W f r is locally maximum The maxima constraints on Wo h x can be decomposed in two conditions 1 At each scale 2 for each local maximum located at x Woih x Wai f 23 2 At each scale 2 the local maxima of Wo h x are located at the abscissa 2 nez Let us first analyze the condition 1 The value of Ws f x at any zo can be written as an inner product in L R Indeed 90 Wa f xo0 f Wei ao J f z Yoi zo x dx thus Was f x0 F Yz o 42 Condition 1 is thus equivalent to Flu vos wi u Alu Yz 2i u Let U be the closure in L R of the space of functions that are linear combinations of functions in the family 43 apoj a onez 44 One can easily prove that the functions that satisfy 43 for all abscissa 2 j n ez are the functions whose orthog onal projection on U is equal to the orthogonal projection of f x on U Let O be the orthogonal complement of U in L R which means that the space O is orthogonal to U and that OSU L R 45 The functions that satisfy 43 for all abscissa 2 j n ez can therefore be written h x f x g x with g x O This defines an affine space that we denote f O If U L R
40. e time it would take to create a singularity at the point xo if we apply a backward heat equation to the signal Let us explain how to measure the smoothing component o as well as the Lipschitz regularity of the underlined singularity We suppose that locally the signal f x is equal to the convolution of a function h a which has a singularity at to with a Gaussian of variance o 1 g xr h g x with g x exp z 32 f z a x Jol Tone Plza 32 We also suppose that h x has a uniform Lipschitz regularity equal to ag in a neighborhood of xo If the wavelet y x is the derivative of a smoothing function z 4 proves that the wavelet transform of f x can be written Wos f t 2 f Bps 0 E hx go 821 2 33 Let us suppose that the function 6 x is close to a Gaussian function in the sense that Oo go x Os x with so V2 07 84 Equation 33 can thus be rewritten a 93 Wi f x 2 h bsa W h x 35 dx So where W h x is the wavelet transform of h x at the scale so Wahl h We 2 This equation proves that the wavelet transform at the scale 2 of a singularity smoothed by a Gaussian of variance g is equal to the wavelet transform of the nonsmoothed singularity h a at the scale so V2 07 Equation 28 of Theorem 1 proves that the Lipschitz regularity is the upper bound of the set of that satisfy Was h x lt K 29 36 This result is valid for
41. e to which we Fig 11 Top left original image Top right image reconstructed from the maxima representation shown in the second column of Fig 9 This reconstruction is performed with FO alternate projections and the SNR is 28 db Lower left image reconstructed from the thresholded modulus maxima shown in the third column of Fig 9 with 10 iterations on the operator P converge is of the same order as the distance between the original image and the image to which we converge Increasing the number of iterations thus does not greatly increase the SNR The top right image in Fig 11 is reconstructed with 1G iterations The SNR is 28 db The reconstructed image has no visual ditference with the original image shown al the top left of Fig 11 which means that the errors are below our visual sensitivity Qualitatively the origmal image is well reproduced because the reconstruction has no spurious oscillation the singularities are not blured and the errors are mostly concentrated at fine scales where our visual sensitivity is nof sO acute The reconstruction algorithm has been tested for a large collection of images including special 2 D functions such as Diracs sinusoidal waves step edges Brownian noises efe For all these experiments the SNR behaves similarly to Fig 10 The visual quality of reconstructed images with 10 iterations is as good as in Fig 11 For image processing applications the numerical precision of this rec
42. eech Signal Processing vol 37 no 12 Dec 1989 12 J Koenderink The structure of images Biol Cybern 1984 13 M Kunt A Ikonomopoulos and M Kocher Second generation image coding technics Proc IEEE vol 74 pp 549 575 Apr 1985 B Logan Information in the zero crossings of band pass signals Bell Syst Tech J vol 56 pp 510 1977 S Mallat Multifrequency channel decompositions of images and wavelet models JEEE Trans Acoustic Speech Signal Processing vol 37 no 12 pp 2091 2110 Dec 1989 Zero crossings of a wavelet transform IEEE Trans Inform Theory vol 37 pp 1019 1033 July 1991 17 S Mallat and W L Hwang Singularity detection and processing with wavelets JEEE Trans Inform Theory vol 38 no 2 Mar 1992 18 D Marr and E Hildreth Theory of edge detection Proc Royal Soc London vol 207 pp 187 217 1980 19 D Marr Vision San Francisco W H Freeman 1982 20 Y Meyer Ondelettes et Operateurs New York Hermann 1990 21 Un contre exemple a la conjecture de Marr et a celle de S Mallat Preprint 1991 22 O Rioul and M Vetterli Wavelets and signal processing IEEE Signal Processing Mag Oct 1991 23 A Rosenfeld and M Thurston Edge and curve detection for visual scene analysis JEEE Trans Comput vol C 20 pp 562 569 1971 24 J Sanz and T Huang Theorem and experiments on ima
43. em is obtained by fixing the parameter y for ej x y and computing the 1 D solution along the x variable between two consecutive modulus maxima The same analysis can be performed on the other component f a y The discrete implementation is thus a straightforward extension of the 1 D algorithm that is applied along the rows and columns of the images that belong to the sequence that we project on I One can verify that if the original image has N pixels the implementation of P p requires O N log N computations In two dimensions we do not introduce any sign constraint as it is done in 1 D reconstructions 117 2 de x y 1 2 2j 243 le z y Dr dx 118 APPENDIX F CONVERGENCE RATE OF THE ALTERNATE PROJECTIONS is a frame of We prove that if v zi x2 z n j EZ U then the alternate projection converges exponentially and we give a lower bound of the convergence rate Let X Wo 9 ez E V and 2 Wo f x Wai g 13 be the error at each modulus maxima location We first prove that there exists a constant C1 gt 0 such that IX PpXl gt Cr e n j EZ 119 731 Let h x jez T Pp Woi9 562 PpX and x hj x W 9 z By definition IX PpX Kele ea gt e jes Past f n j EZ n B Pa ts 120 We saw in Appendix E that e x satisfies the differential equation 112 therefore by integrating by parts we obtain
44. ersity New York NY 10012 He is now with the Department of Mathematics University of California Los Angeles CA IEEE Log Number 9200110 wavelet transform maxima across scales Lipschitz exponents and smoothing factors are numerical descriptors that allow us to discriminate the intensity profiles of different types of edges An important open problem in computer vision is to under stand how much information is carried by multiscale edges and how stable a multiscale edge representation is This issue is important in pattern recognition where one needs to know whether some interesting information is lost when representing a pattern with edges We study the reconstruction of 1 D and 2 D signals from multiscale edges detected by the wavelet transform modulus maxima It has been conjectured 16 18 that multiscale edges characterize uniquely 1 D and 2 D signals but recently Meyer 21 has found counterexamples to these conjectures In spite of these counterexamples we show that one can reconstruct a close approximation of the original signal from multiscale edges The reconstruction algorithm is based on alternate projections We prove its convergence and derive a lower bound for the convergence rate Numerical results are given both for 1 D and 2 D signals The differences between the original and reconstructed images are not visible on a high quality video monitor The ability to reconstruct images from multiscale edges has many applic
45. es slightly at the location where Wz f x has this sharp change These oscillations are similar to a Gibbs phenomenon Appendix E explains how to modify the alternate projections in order to suppress these oscillations Numerical experiments show that this oscillation removal does not perturbate the convergence of the algorithm C Numerical Reconstruction of 1 D Signals from Local Maxima There are several open issues behind the reconstruction algorithm that we described From the results of Meyer s work 21 we know that in general we cannot reconstruct exactly a function from the modulus maxima of its wavelet transform Our algorithm approximates this inverse problem by replacing the maxima constraint by the minimization of a norm that yields a unique solution We thus do not converge 719 Fig 4 Approximation of the wavelet transform of f x is reconstructed by alternating orthogonal projections on an affine space T and on the space V of all dyadic wavelet transforms The projections begin from the zero element and converge to its orthogonal projection on TN V toward the wavelet transform of the original signal but toward some other wavelet transform that we hope to be close to the original one We also explained that the computation of the solution might be unstable in which case the alternate projections converge very slowly It is therefore important to measure how far we are from the convergence point after a given number
46. f each of these modulus maxima as well the values of the modulus My f x y and the angle As f x y at the corresponding locations The circle image at the top of Fig 8 has 128 by 128 pixels The first two columns of Fig 8 give the discrete wavelet transform W f and we f for 1 lt j lt 8 The next section explains how to define such a discrete dyadic wavelet transform and how to solve border problems The reader that is not interested by numerical implementations can skip Section VI B The discrete modulus images M4 f and angle images AS f are shown along the next two columns Along the border of the circle the angle value turns from 0 to 27 and the modulus has a maximum amplitude When the scale 2 is larger than 2 we see that the circle is deformed due to the image periodization that we use for border computations The position of the modulus maxima at all scales is given in the last column on the right The original Lena image is shown at the top left of Fig 12 and has 256 by 256 pixels The first column of Fig 9 displays its discrete modulus images M3 f and the second column gives the position of the modulus maxima for 1 lt j lt 9 At fine scales there are many maxima created by the image noise At these locations the modulus value W f x y W2 f x y 67 argument W3 f x y iWS f x y 68 ton TREE ks Fig 8 from the left show respectively W Original image at the lop ns oy y 128 pi
47. finer scale available and the Lipschitz regularity is computed by finding the coefficient o such that K 2 approximates at best the decay of W f x over a given range of scales larger than 1 In Fig 3 b in the neighborhood of x 2 the maxima values of W f x remain constant over a large range of scales Equation 31 implies that the Lipschitz regularity ap is equal 715 to 0 at that point which means that this singularity is a discontinuity In the edge detection procedure described in Section II we only keep the local maxima of the wavelet transform modulus It has been proved 17 that if a signal is singular at a point xo there exists a sequence of wavelet transform modulus maxima that converge to xo when the scale decreases Hence we detect all the singularities from the positions of the wavelet transform modulus maxima Moreover the decay of the wavelet transform is bounded by the decay of these modulus maxima and we can thus measure the local uniform Lipschitz regularity from this decay A signal is often not singular in the neighborhood of local sharp variations An example is the smooth edge at the abscissa 1 of Fig 3 a It is generally important to estimate the smoothness of the signal variation in such cases We model a smooth variation at ro as a singularity convolved with a Gaussian of variance g Since the Gaussian is the Green s function of the heat equation one can prove that g is proportional to th
48. ge reconstruc tion from zero crossings Res Rep RJ5460 IBM Jan 1987 25 A Witkin Scale space filtering in Proc Int Joint Conf Artificial Intell 1983 D Youla and H Webb Image restoration by the method of convex projections IEEE Trans Med Imaging vol MI 1 pp 81 101 Oct 1982 27 A Yuille and T Poggio Scaling theorems for zero crossings IEEE Trans Patt Anal Machine Intell vol PAMI 8 Jan 1986 28 Y Zeevi and D Rotem Image reconstruction from zero crossings IEEE Acoustic Speech Signal Processing vol 34 pp 1269 1277 1986 29 S Zhong Edges representation from wavelet transform maxima Ph D thesis New York Univ Sept 1990 10 11 14 15 16 26 Stephane Mallat was born in Paris France He graduated from Ecole Polytechnique in 1984 and from Ecole Superieure des Telecommunications Paris in 1985 He received the Ph D degree in electrical engineering from the University of Pennsylvania Philadelphia in 1988 Since September 1988 he has been Assistant Professor in the Computer Science Department of the Courant Institute of Mathematical Sciences New York University New York NY His research interests include computer vision signal processing and applied mathematics Dr Mallat received the 1990 IEEE Signal Processing Society s paper award Sifen Zhong received the B A degree in 1985 in mathematics and computer sc
49. high frequencies The numerical precision of reconstructions 1s thus not improved with this other wavelet A discrete analysis of the completeness conjecture was done independently by Berman 1 who found numerical examples that contradict the completeness conjecture We explained in Section II that for a wavelet equal to the first derivative of a smoothing function the local minima of the 717 wavelet transform modulus correspond to slow variation points of the signal Hence among all the wavelet transform extrema we detect only the points where the wavelet transform modulus is locally maximum For the quadratic wavelet of Fig 1 a since the wavelet transform local extrema do not provide a complete signal representation the subset of modulus maxima is certainly not complete either The next section describes an algorithm that still recovers a precise approximation of the original signals from these modulus maxima B Reconstruction Algorithm Let f x L R and Wo f x ez be its dyadic wavelet transform We describe an algorithm that reconstructs an approximation of Wa f x z given the positions of the local maxima of W3 f x and the values of Ws f a at these locations For this purpose we characterize the set of functions h x such that at each scale 27 the modulus maxima of Wo h x are the same as the modulus maxima of o f x We suppose that the wavelet y x is differentiable in the sense of Sobolev Sinc
50. ience from Queen s College New York NY and the M S and Ph D de grees in computer science from the Courant Institute of Mathematical Sciences New York University New York NY in 1987 and 1991 respectively He worked for two years as a research assistant at the New York University Robotics Research Labo ratory Presently he is a visiting assistant researcher at the Department of Mathematics of the University of California Los Angeles He will join Cognitech Inc an image processing company in Santa Monica CA in 1992
51. ignal smoothed at the scale 2 Fig 2 c gives the locations and values of the local maxima of the dyadic wavelet transform modulus as in a Canny edge detection At each scale 27 each modulus maximum is represented by a Dirac that has the same location and whose amplitude is equal to the value of Wz f x The modulus maxima detection is an adaptive sampling that finds the signal sharp variation points B Discrete Wavelet Transform In numerical applications the input signal is measured at a finite resolution therefore we cannot compute the wavelet transform at an arbitrary fine scale Let us normalize the finest scale to 1 In order to model this scale limitation we introduce a real function 2 whose Fourier transform is an aggregation of 7 24w and 2 w at scales 27 larger than 1 90 elo So P 2 w X 2 w 20 j l We suppose here that the reconstructing wavelet x w is such that 1 w w is a positive real even function One can prove that property 13 implies that the integral of x is equal to 1 and hence that it is a smoothing function Let Sj be the smoothing operator defined by T S fla f o 2 with do 2 AS CD If the scale 27 is larger the more details of f x are removed by S For any scale 27 gt 1 20 yields 713 Sif 92 A ee SE SEER ee as ee ee eee ee Nee 93 r E Ce ane ee EY Vee CTY Lee ae eT eee eee oe a AA ye gt 2 A gt _ _1____ 9 _1t _ ___
52. ive A tempered distribution f x has a uniform Lipschitz regularity equal to o over Ja b if and only if for any a lt ag there exists K such that 30 Wz f x lt K 2 31 Since the Lipschitz regularity of a Dirac is 1 this result implies that the maxima values of W 5 x increase propor tionally to the scale 2 This can indeed be verified in Fig 3 b In practice we can only process discrete signals that ap proximate the original function at a finite resolution which we normalize to 1 Strictly speaking it is not meaningful to speak about singularities discontinuities or Diracs In fact we cannot compute the wavelet transform at scales finer than 1 and thus cannot verify 31 at scales smaller than 1 Even though we are limited by the resolution of measurements we can still use the mathematical tools that differentiate singularities Suppose that the approximation of f z at the resolution 1 is given by a set of samples fn ez with fn 0 for n lt no and f 1 for n gt no like at abscissa 2 of Fig 3 a At resolution 1 f x behaves as if it has a discontinuity at n no although f x might be continuous at no with a continuous sharp transition at that point which is not visible at this resolution The characterization of singularities from the decay of the wavelet transform gives a precise meaning to this discontinuity at the resolution 1 We measure the decay of the wavelet transform up to the
53. ize this sum we minimize separately each component de leal 2 1 eH Let zo and z be the abscissa of two consecutive modulus maxima of W fz x Since h 7 T we have 20 Wai f 0 9 x0 e 21 Was f w1 g 21 Between the abscissa xq and x the minimization of 109 is equivalent to the minimization of 7 de x cor 21 ae 110 111 The Euler equation associated with this minimization is d x 2 j e z 2 j dr2 for x x o x The constraints 110 are the border conditions of this membrane equation The solution is 0 112 g r ae be 113 where the constants a and 2 are adjusted to satisfy equations 110 In numerical computations Wwe f is a uniform sampling of Wo f x at the rate 1 and has a total of N samples At cach scale 2 the operator P p modifies a ares signal gf 95 1 lt n lt n by adding a discrete signal ef ej n i lt n lt N that is computed from 113 between two consecutive modulus maxima This requires O N computations Since there are at most log N 1 scales the total number of computations to implement Pp is O N logs N We know that the modulus maxima of the original wavelet transform are only located at the positions z We can thus also impose sign constraints in order to suppress any spurious oscillation in the reconstructed wavelet transform This is done by imposing that the solution be
54. l signal as in Fig 6 b the error is measured with respect to the signal to which we converge After 30 iterations the slope of the SNR curve is constant which proves that the convergence is exponential but the convergence rate is slow In Fig 6 b and Fig 7 the increase of the SNR slows down after approximately 20 iterations At this point the distance between the reconstructed signal and the signal to which we converge is of the same order as the distance between the original signal and the signal to which we converge Increasing the number of iterations slowly reduces the distance with respect to the point to which we converge but does not largely decrease the distance with respect to the original signal This is why the SNR in Fig 7 continues to increase slowly whereas the SNR in Fig 6 b reaches a maximum on the order of 38 db We made extensive numerical tests including reconstruc tions of special functions such as sinusoidal waves Gaussians step edges Diracs fractals and the counter example of Meyer given by 41 In all these examples the SNR has the same type of behavior as in Figs 6 and 7 In most cases after 30 iterations the relative increase of precision that is obtained 60 i 50 2 40 7 93 30 92 91 20 10 0 20 40 60 80 100 a 40 35 30 25 20 15 10 0 20 40 60 80 100 Fig 6 a SNR for the reconstruction of the wavelet transform we fasa function of the number of iterations on the operator P
55. lem in image processing is to code images with a minimum number of bits for transmission or storage To obtain high compression rates in image coding we cannot afford to code all the information available in the image It is necessary to remove the parts of the image components that are not important for visualization A major problem is to identify the important information that we need to keep From this point of view the problems encountered in compact image coding are similar to computer vision tasks where one also wants to extract the important information for recognition purposes Since edges provide meaningful features for image interpretation it is natural to represent the image information with an edge based representation in order to select the information to be coded Previous edge coding algorithms have already been developed by Carlsson 3 and Kunt et al 13 but at a single scale This section describes a compact coding algorithm based on the wavelet transform modulus maxima The coding algorithm involves two steps First we select the edge points that we consider important for the visual image quality This preprocessing is identical to the feature extraction stage of a pattern recognition algorithm We then make an efficient coding of this edge information Selection of the most important edge curves can require sophisticated algorithms if we take into account the image context For example in the Lena image
56. longs to an appropriate convex set Y Let sign x be the sign of the real number z Let Y be the set of sequences g x jez K such that for any pair of consecutive maxima positions x3 7 a and x z2 27 4 sign g z sign z3 if sign x3 sign x 1 49 2 sign sign ry 41 z4 if sign x sign x 41 The set Y is a closed convex and Wz f ez Y Instead of minimizing over T N V as explained in Section V B we can minimize it over Y NI N V We thus alternate projections on Y F and V To compute the orthogonal projection on the convex Y we need to solve an elastic membrane problem under constraints This can be done with an iterative algorithm that is computationally intensive Instead we implement a simpler projector Py on Y which is not orthogonal with respect to the norm Let g jez E K and Py 9 x lt 7 hj z For each index j h x is obtained by clipping the oscillations of g z If the original signal has N samples at each scale 2 the discrete implemen tation of the clipping procedure requires O N computations There are at most log N 1 scales therefore the total number of computations to implement Py is O N log N Since this projector Py is not orthogonal the iteration on the alternate projections operator P PyoPpoPy is not guaranteed to converge Numerical experiments show that in most cases after a few iterations we stay inside Y even
57. ly modulus maxima of the wavelet transform but that they minimize a Sobolev norm defined by 2 pr Wo A x y W3 A z y ez 00 uaa wn j 00 WLR aW2 h 2j PTAA gie no z 2 EA I 85 The minimization of this norm creates a wavelet transform whose horizontal and vertical components have an L R norm that is as small as possible In conjunction with condition 1 this has a tendency to create modulus maxima at the po sitions z y The partial derivative components are added in order to create a wavelet transform with as few spurious oscillations as possible Since W h x y is computed by smoothing the signal and taking the partial derivative along zx it oscillates mostly along the x direction and we use a partial derivative along x in 85 to minimize these oscillations The transpose result is valid for W2 h x y The weight on the derivative components is proportional to the scale 2 because the smoothness W4 h x y and W2 h a y increases with the scale 2 Let y3 z y Ov zw and yt z y ee If there exist two constants A gt 0 and Bs gt O such that for all ws wy R o Ass Y Pwe Huy we Puy j MALLAT AND ZHONG CHARACTERIZATION OF SIGNALS 00 D Pwe uy we Mwy P lt Bs J 86 then for any function h z y L R the norm defined in 85 is equivalent to the L R norm AsllAll lt IAI lt Bs lAll 87 Similar to 49
58. m the scale 2 up to the scale 2 and converge to the abscissa xo Let a be the value of the wavelet transform at the maximum location at the scale 27 and let us also suppose that in a given neighborhood of xo the wavelet transform modulus is smaller than a This means that the signal change at xo is the sharpest variation in this neighborhood We compute the three values K o and ap so that the inequality of 37 is as close as possible to an equality for each maximum a These values are obtained by minimizing a 1 2 loga o 2 38 I X 108 la loga K j j 1 This is done with a steepest gradient descent algorithm The value K gives the amplitude of the sharp variation When computing the values of o and a from the evolution of the maxima across scales in Fig 3 we have a numerical error of less than 10 which is mainly due to the fact that the wavelet we use is not the derivative of a Gaussian but is only an approximation In this case 0 x is the cubic spline shown in Fig 1 b When the variance o increases the measurement of o becomes more unstable because the smoothing removes the fine scale components that characterize reliably ag For singularities of fractal textures such as in the right part of Fig 2 a this analysis is not valid because singularities are not isolated and none of the singularities dominate the others in a given neighborhood The behavior of the wavelet transform modulus
59. n Section V C A 2 D extension of the wavelet transform is given in Section VI A 0162 8828 92 03 00 1992 IEEE MALLAT AND ZHONG CHARACTERIZATION OF SIGNALS and its discrete version is explained in Section VI B Fast 2 D wavelet algorithms are given in Appendix D Section VII differentiates the edges of an image from the evolution across scales of the wavelet modulus maxima The reconstruction of images from multiscale edges is explained in Section VIII A and numerical examples are shown in Section VIII B Section IX describes an application to compact image coding Notation L R denotes the Hilbert space of measurable square integrable 1 D functions f x For f L R and g L R the inner product of f x with g x is written 00 g x f Jai co The norm of f x L R is given by 2 E 2 I f Ode We denote the convolution of two functions f x L R and g x L R by 00 f g 2 J Fee eee oO The Fourier transform of f x L R is written f w and is defined by 20 fw f Haje tda L R is the Hilbert space of measurable square integrable 2 D functions f x y The norm of f x y L R is given by 2 a re 2 IfI l J Fe dedy The Fourier transform of f x y L R is written A f wz wy and is defined by 4 00 00 f we Wy j f f a yje tte dady IJ MULTISCALE EDGE DETECTION Most multiscale edge detectors smooth the signal at various
60. nd only if there exist two constants A3 gt 0 and B3 such that for any function g U Asla lt S gt 2 g u v2 a4 u lt Ballgll GL n JEZ When the two constants Ay and B of 48 are different the norm is not equal but is equivalent to the classical L R norm In this case the stability also depends on whether the family of wavelets is a frame of U The closer to 0 the value of aoa the more stable the computations Outside of a few particular cases it is difficult to prove analytically whether a given family of wavelets v Palri x INE not a frame of the space U that it generates because the points xi are not uniformly distributed Let us now describe an algorithm that computes the solu tion of our minimization problem Instead of computing the solution itself we reconstruct its wavelet transform with an algorithm based on alternate projections Let K be the space of all sequences of functions io jez such that iS or is Z2 o 2 j 2 e Do igi e lt 400 62 The norm defines a Hilbert structure over K Let V be the space of all dyadic wavelet transforms of functions in L R Equation 49 proves that V is included in K Let I be the affine space of sequences of functions g x z K such that for any index j and all maxima positions x7 ga Wy F One can prove that I is closed in K The dyadic wavelet trans forms that satisfy condition 1 are the sequence
61. oarse scale 27 the sequence of discrete signals Sirf Ws ie is called the discrete dyadic wavelet transform of D Sif n 7 The coarse signal S4 provides the signal com ponents below the scale 27 A fast discrete wavelet transform algorithm and its inverse are described in Appendix B In practice the original discrete signal D has a finite number N of nonzero values D dn ic lt n To solve the border problems we use the same periodization technique as in a cosine transform We suppose that our signal has a period of 2N samples and we extend it with a symmetry for N lt n lt 2N dn don41 n By periodizing the signal with a symmetry we avoid creating a discontinuity at the borders The discrete wavelet coefficients are also 2N periodic If the wavelet is antisymmetrical with respect to 0 as in Fig 1 a the wavelet coefficients are antisymmetrical at the borders For the class of wavelets defined in Appendix A one can also prove that when the scale is as large as the period 2 2N S f is constant and equal to the mean value of the original signal D We thus decompose any signal of N samples over J loga N 1 scales Appendix B describes a fast discrete wavelet transform algorithm that requires O N log V operations The fast inverse wavelet transform also requires O N log V operations From the discrete wavelet transform at each scale 25 we detect the modulus maxima by finding the points where W f n
62. odulus maxima are the points of the modulus images M f that are larger than the two neighbors whose positions are in the direction indicated by the angle value of Ag f at the corresponding location VII CHARACTERIZATION OF IMAGE EDGES Sharp variations of 2 D signals are often not isolated but belong to curves in the image plane Along these curves the image intensity can be singular in one direction while varying smoothly in the perpendicular direction It is well known that such curves are more meaningful than edge points by themselves because they generally are the boundaries of the image structures For discrete images we reorganize the maxima representation into chains of local maxima to recover these edge curves As in one dimension we then characterize the properties of edges from the modulus maxima evolution across scales At a scale 27 the wavelet modulus maxima detect the sharp variation points of f 9 x y Some of these modulus maxima define smooth curves in the image plane along which the profile of the image intensity varies smoothly At any point along a maxima curve V f 49 z y is perpendicular to the tangent of the edge curve We thus chain two adjacent local maxima if their respective position is perpendicular to the direction indicated by As f x y Since we want to recover edge curves along which the image profile varies smoothly we only chain together maxima points where the modulus Mo f x y has close v
63. onstruction algorithm is sufficient even if we limit the number of iterations below 10 Since each iteration requires O N log N computations this reconstruction can be implemented in hardware for real time applications The reconstruction algorithm is stable for precisions on the order of 30 db We can therefore slightly per turbate the wavelet transform modulus maxima and reconstruct a close image The lower left image in Fig 11 is reconstructed from the modulus maxima shown in the third column of Fig 9 By thresholding the wavelet transform modulus maxima based on their modulus values we suppressed the modulus maxima produced by the image noise and the light textures As expected these textures have disappeared in the reconstructed image but the sharp variations are not affected In the lady s shoulder the thresholding removes the maxima created by Patt IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE VOL 14 NO 7 HILY 1882 the image noise and the reconstructed image reproduces a skin quality that is much smoother whereas the boundaries of the shoulder are kept sharp This thresholding can be viewed as a nonlinear noise removal technique Hwang and Mallat 17 have developed a more sophisticated procedure to suppress white noises from images which removes the maxima produced by the noise through an analysis of their behavior across scales EX Compact IMAGE CODING FROM MULTISCALE EDGES An important prob
64. res are generally much less visible than distortions of edges and a separate coding of these two types of features can be adapted to the specificity of the visual perception X CONCLUSION We showed that multiscale edges can be detected and characterized fram the local maxima of a wavelet transform a Fig 12 Top left original image of 246 by 256 pixels Top right recon structed image from the coded multiscale edge representation Image 4 requires O 30 b per pixel and image and b requires 0 24 b per pixel Bottom left image f subsampled by a factor 2 along its rows and columss Bottom right position of the modulus maxima selected and encoded at the scales 2 2 ee ap 93 Fe YEG ay ms One can estimate the Lipschitz regularity as well as the smoothing component of sharp variation points from the evolution of the wavelet maxima across scales We believe that this complement of information is important for pattern recognition algorithms based on edges The reconstruction algorithms that are described in 1 D and 2 D recover a close approximation of the original signals For images the reconstruction errors are below our visual sensitivity and can thus be neglected in image processing or computer vision applications To reconstruct such signals requires few iterations that can be implemented in real time on a pipeline hardware architecture As an application we 728 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE
65. s of functions that belong to A VOL We must therefore find the element of A whose norm is minimum This is done by alternating projections on V and I Equation 17 shows that any dyadic wavelet transform is invariant under the operator Py Wow 53 For any sequence X g 7 K it is clear that Py X V therefore Py is a projector on V We saw in 19 that this operator is characterized by the kernels K j xo Wo x One can easily prove that the projector Py is self adjoint and therefore orthogonal if and only if the kernels Ay x are symmetrical functions This is the case if the wavelet x is symmetrical or antisymmetrical For the wavelet shown in Fig l a the orthogonal projection on the space V is thus implemented by applying the operator W followed by the operator W The fast discrete implementation of these operators is given in Appendix B Appendix E characterizes MALLAT AND ZHONG CHARACTERIZATION OF SIGNALS the projection on the affine set I which is orthogonal with respect to the norm We prove that this operator Pr is implemented by adding piecewise exponential curves to each function of the sequence that we project on I Let P PyoPp be alternate projections on both spaces Let P be n iterations over the operator P Since T is an affine space and V a Hilbert space a classical result on alternate projections 26 proves that for any sequence of functions X 9i 2 je
66. singularity at the abscissa 3 of Fig 3 b produces wavelet 29 MALLAT AND ZHONG CHARACTERIZATION OF SIGNALS transform maxima that increase when the scale decreases Such singularities can be described with negative Lipschitz exponents which means that they are more singular than discontinuities The signal is then viewed as a tempered distribution At the abscissa 3 of Fig 3 b this distribution is locally equal to a Dirac The reader might want to consult Folland 7 for a quick presentation of the mathematical theory of distributions The wavelet transform of tempered distributions is well defined if the wavelet y x is smooth enough 17 For example if y x is continuous the wavelet transform of a Dirac 6 z is given by Was 6 x 6 Woi x Ya x To extend Lipschitz exponents to distributions we say that a distribution has a uniform Lipschitz regularity equal to a on ja b if and only if its primitive has a uniform Lipschitz regularity equal to a 1 on Ja b For example the primitive of a Dirac centered at xo is a function that is bounded and has a discontinuity at zo step edge The uniform Lipschitz regularity of the primitive of this Dirac is thus equal to 0 in the neighborhood of xo Hence a Dirac centered at xo has a uniformly Lipschitz regularity equal to 1 in the neighborhood of zo One can prove that Theorem 1 is also valid for negative Lipschitz exponents Let ag lt 1 be a real number that may be negat
67. ss ca mre EEO A gi 2 wyk Daw d wy 69 We shall suppose that x y is real From the admissibility condition 61 one can derive that the integral of Ur y is equal to 1 which means that it is a smoothing function We define the smoothing operator S by 5 ane l ee ee A koft WA So fira foaie y with oxir y ef e DN i H 4 i I ee oe yi ith 3 A i uy pa Aj i pe 70 The the scales 1 and 2 lt y provides the details al ar wavelet transform between yh e Wh flay W2 fla y Se available in S fir y but a have disappeared in Sos f x y TRANSACTIONS ON PATTERN ANALYSES AND MAC HINE INTELLIGENCE VOL 4 NOL 7 JULY 1992 Fig 9 Original Lena image is at the top left of Fig 1 and has 256 by 286 pixels The first column gives Hs modulus images aT i o The second column displays the position of the local maxima of Afs f i lt j lt 9 The third column gives the positions of local maxima where the modulus value is larger than a given threshold Local maxima that correspond to Hight texture variations are removed by the thresholding At the output of a camera digitizer an ee is a finite energy 2 D discrete signal D dion inmez Similar to the 1 D case one can prove that here exists a Function f z e DIR not ie such that af reer Vo Vn The discrete image ean thus be rewritten Prt ijin mez For the class of wavelets ytGr
68. t caci scale 27 the algorithm decomposes s f into Sat a wat f and Woe j 0 while j lt J Word 94 f Gj D tid E no Ne Gj So aad m S5 f Hj H j j 1 end of while The proof of this algorithm is based on the properties of the wavelets 71 z y and y x y described in Appendix C If the original image S1 f n m _mycz2 has N 2 nonzero pixels the complexity of the algorithm is O N log N As explained in Section VI B border problems are solved by making a symmetry of the image with respect to each of its borders and a periodization The separable convolutions must take into account this border procedure As in the 1 D case the reconstruction algorithm computes S f by reconstructing at each scale 2 the signal Sos W from S f W4 df and W2 df The complexity of this feconsiruction algorithm is also O N log N f a while i j gt 0 daf Ay Wa S Kya Lj Ay Woy f Lj A 1 Sd f Fi 1 H 1 J 9 1 end of while APPENDIX E PROJECTION OPERATOR ON F In this appendix we characterize the orthogonal projection on I in one and two dimensions and explain how to suppress oscillations for 1 D reconstructions We first study the 1 D case The operator P p transforms any sequence 9 2 7 K into the closest sequence h x T with respect to the norm Let e x h x g x Each function h z is chosen so that de 3 lejl 2 l a 109 j oo is minimum To minim
69. two strictly positive constants A and By such that 00 WweRAr lt X h 2 w lt B 12 j we ensure that the whole frequency axis is covered by dilations of w by 2e so that f w and thus f x can be re covered from its adie wavelet transform The reconstructing wavelet x x is any function whose Fourier transform satisfies 00 S o 24w X 24w 1 j 00 13 If property 12 is valid there exists an infinite number of func tions x w that satisfy 13 The function f x is recovered from its dyadic wavelet transform with the summation 00 S0 Wai f x2 2 j 00 14 This equation is proved by computing its Fourier transform and inserting 11 and 13 With the Parseval theorem we derive from 11 and 12 a norm equivalence relation 00 gt Was fe lt Bill fl j Aillf l lt 15 This proves that the dyadic wavelet transform is not only complete but stable as well If Ei is closer to 1 it will be more stable A dyadic wavelet transform is more than complete it is redundant Any sequence g lt z with g x L R is not necessarily the dyadic wavelet transform of some function in L R We denote by W the operator defined by gt gj X2i Z j wt 9 x jez 16 The reconstruction formula 14 shows that gi x T jez is the dyadic wavelet transform of some function in L R if and only if W W 17 g t teZ 9 2 sez If we
70. urves have been called finger prints by Witkin 25 We say that a modulus maxima at the scale 2 propagates to a maxima at the coarser scale 2 1 if and only if both maxima belong to the same maxima curve in the scale space plane s r In Fig 3 there is one sequence of maxima that belongs to the same maxima curve and converges to the position of the discontinuity at x 2 For the Dirac at abscissa 3 there are two such sequences Each one gives information respectively on the left and the right part of the Dirac singularity In order to find which maxima propagate to the next scale one should compute the wavelet transform on a dense sequence of scales However with a simple ad hoc algorithm one can still estimate which maxima propagate to the next scale by looking at their value and position with respect to other maxima at the next scale The propagation algorithm supposes that a modulus maximum propagates from a scale 2 to a coarser scale 2 1 if it has a large amplitude and if its position is close to a maximum at the scale 27 1 that has the same sign This algorithm is not exact but saves computations since we do not need to compute the wavelet transform at any other scale The Lipschitz regularity as well as the smoothing variance o of a sharp variation point are then computed from the evolution of the modulus maxima that propagate across scales Let us suppose that we have a sequence of modulus maxima that propagate fro
71. uter vision have introduced the concept of multiscale edge detection 18 23 25 The scale defines the size of the neighborhood where the signal changes are computed The wavelet transform is closely related to multiscale edge detection and can provide a deeper understanding of these algorithms We concentrate on the Canny edge detector 2 which is equivalent to finding the local maxima of a wavelet transform modulus There are many different types of sharp variation points in images Edges created by occlusions shadows highlights roofs textures etc have very different local intensity profiles To label more precisely an edge that has been detected it is necessary to analyze its local properties In mathematics singularities are generally characterized by their Lipschitz exponents The wavelet theory proves that these Lipschitz exponents can be computed from the evolution across scales of the wavelet transform modulus maxima We derive a numerical procedure to measure these exponents If an edge is smooth we can also estimate how smooth it is from the decay of the Manuscript received January 10 1990 revised January 29 1992 This work was supported by the NSF grant IRI 890331 AFOSR grant AFOSR 90 0040 and ONR grant N00014 91 J 1967 Recommended for acceptance by Associate Editor R Woodham S Mallat is with the Courant Institute New York University New York NY 10012 S Zhong was with the Courant Institute New York Univ
72. vative of the signal smoothed at the scale s The local extrema of W f z thus correspond to the zero crossings of W f x and to the inflection points of f a In the particular case where x is a Gaussian the zero crossing detection is equivalent to a Marr Hildreth 19 edge detection whereas the extrema detection corresponds to a Canny 2 edge detection When the scale s is large the convolution with removes small signal fluctuations we therefore only detect the sharp variations of large structures Detecting zero crossings or local extrema are similar pro cedures but the local extrema approach has some important advantages An inflection point of f x can either be a maximum or a minimum of the absolute value of its first derivative The maxima of the absolute value of the first derivative are sharp variation points of f s x whereas the minima correspond to slow variations With a second derivative operator it is difficult to distinguish these two types of zero crossings On the contrary with a first order derivative we easily select the sharp variation points by detecting only the local maxima of W f x In addition zero crossings give position information but do not differentiate small amplitude fluctuations from important discontinuities When detecting local maxima we can also record the values of W f x at the maxima locations which measure the derivative at the inflection points Section IV explains
73. w is larger than its two closest neighbor values and strictly larger than at least one of them We record the abscissa n w and the value W2 f n w at the corresponding locations 26 IV ANALYSIS OF THE MULTISCALE INFORMATION One signal sharp variation produces modulus maxima at different scales 27 We know that the value of a modulus maximum at a scale 2 measures the derivative of the signal smoothed at the scale 2 but it is not clear how to combine these different values to characterize the signal variation The wavelet theory gives an answer to this question by showing that the evolution across scales of the wavelet transform depends on the local Lipschitz regularity of the signal This section explains what a Lipschitz exponent is and how this exponent is computed from the wavelet transform maxima S4 f Sn 0 1 2 a 3 4 5 Waf Wf Waf A Waf N b Fig 3 a Four sharp variation points of this signal have a different Lipschitz regularity ao and smoothing variance g These values are given respectively by ao 0 0 3 ag 0 0 0 ao l o 0 and ao l g 4 b behavior of the modulus maxima across scales depends on the Lipschitz regularity ag and the smoothing factor A more detailed mathematical and numerical analysis of this topic can be found in 17 When the signal is not singular we show that one can still measure how smooth the signal is by estimating the decay of the w
74. ws the increase of the SNR which is computed with respect to the original signal This SNR is an aggregation of the wavelet transform SNR at all scales The signal in Fig 5 b is reconstructed by applying the inverse wavelet operator on the reconstructed wavelet transform after 20 iterations In this case the SNR is 34 6 db The remaining error after n iterations has two components The first one is the distance between the reconstructed wavelet transform and the wavelet transform to which we converge The other one is the distance between the wavelet transform to which we converge and the wavelet transform of the original signal We saw that the convergence rate of the algorithm is related to the frame properties of the family of wavelets 720 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE VOL 14 NO 7 JULY 1992 a b Fig 5 a Original signal b signal reconstructed with 20 iterations from the modulus maxima shown in Fig 2 c defined by the maxima positions In numerical computations there is a finite number of maxima therefore the family of wavelets that generates U is finite A finite family of vectors is always a frame but the frame bound A can be very smali The lower bound of the convergence rate given by 55 can thus also be very small Fig 7 is the SNR of the reconstructed signal computed with respect to the signal to which we converge Instead of measuring the error with respect to the origina
75. x which are defined with respect to the wavelet 7 2 ae From 4 5 and 39 we derive that Cn Wait 4a 1 Wo fa 40 To record the zero crossing positions and integral values of W2 f x is therefore equivalent to recording the positions where W3 f x has local extrema and the value of W f x at the corresponding locations Meyer 21 proved that the completeness of this representation depends on the choice of the smoothing function x but that the conjecture is not valid in general Let fo x i 1 cos x if z Ea 41 0 otherwise For the wavelet shown in Fig 1 a Meyer 21 found a noncountable family of functions Fez fo x pa Xe z such that at all scales 27 Wo fe x and Wo fo x have the same extrema positions and values The functions x are small high frequency perturbations which are implicitly defined by constraint equations that guarantee that the local extrema of W3 fo x are not modified It seems that in order to maintain the local extrema of Wo fo x unchanged the perturbations xe x must remain small which would explain the quality of the signal reconstructions obtained in 16 but this has not been proved For another wavelet defined by w r ate with x folz Meyer proved that any function of compact support is uniquely characterized by the zero crossings and integral values of its dyadic wavelet transform This characterization is however not stable at
76. xamples of images coded with this algorithm The same length and amplitude thresholds were used for each of these images to select the edge chains at the scale 2 For each example we display at the top left the original image at the top right the image reconstructed from the coded representation at the bottom right the edge map at the scale 27 that is encoded and at the bottom left the subsampling of the low frequency image Se f Each original image has 256 by 256 pixels The total amount of data to code the reconstructed images is 0 30 b per pixel for Fig 12 a and 0 24 b per pixel for Fig 12 b The compression rate varies with the number of edge points that remain after the selection operation This type of coding removes the image textures however H does not produce distortions such as Gibbs phenomena For the Lena image errors are particularly visible around the eyes because too many edge points have been removed in this region by our simple selection algorithm Although a lot of details have been removed in the coded images they remain sharp and most of the information is kept This compact coding algorithm is a feasibility study and if can certainly be improved both at the selection and the coding stages Por applications to images where textures are important Froment and Mallat 8 have extended the method by developing an algorithm that makes a specifie coding of textures after this edge based coding Distortions of textu
77. xels The firsi two columa s and oi p Ehe lt a i 2 p lt PR scale increases from top to bottom Blac fe grey and white pixels indicate respectively negative zero and positive pixel values The third column 7 rh arf ei j I lt values Whereas white ones coe to Highest values The fourth column displays the modulus images Black pixels indicate zero shows the angle images 4 48 F The angle value turns from 0 Ow hite to 2a black along the citcle contour The images of the last column display in white the pomts where i Be 3 f has local maxima in the direction indicated wed has a small amplitude The third column displays the maxima whose modulus is larger than a given threshold at all scales The edge points with a high modulus value correspond to the sharper intensity variations of the image At coarse scales the modulus maxima have different positions than at fine scales This is due to the smoothing of the image by fy 2 y B Discrete Wavelet Transform of Images Images are measured at a finite resolution therefore we cannot compute the wavelet transform at scales below the limit set by this resolution As in one dimension in order to model the limitation of resolution we introduce a smoothing function x y whose Fourier transform is an aggregation of the wavelet components dilated by scales larger than wie KT F iZ a i 4 yni j low wey My POO ab Dor 2 wy JX Hwa aly pis
78. z ek lim POX P 4X n o 54 Alternate projections on and V converge strongly to the orthogonal projection on A If X is the zero element of K which means that g x O for all 7 Z the alternate projections converge to the element of A which is the closest to zero and thus whose norm is minimum This is illustrated by Fig 4 This iterative algorithm can be related to techniques based on frame operators for reconstructing signals from irregular samplings 9 If the minimization problem is unstable the conver gence of the alternate projections is extremely slow We saw that the numerical stability depends on whether Vas x2 2 is a frame of U Appendix F ZN EA proves that if Vha xi eg is a frame of U n jJ EL and if there exists a constant 0 lt D lt 1 such that at all scales 2 the distances between any two consecutive maxima satisfy x _ Taai gt D2 then the convergence is exponential Moreover there exists a constant R such that for any X K Pp P lt n 2 X P4X lt RO 52 DAs 55 2 where the A is the frame bound defined in 51 and B the norm equivalence bound defined in 48 This equation gives a lower bound for the convergence rate and shows how it decreases when the frame bound Ag goes to zero When the original wavelet transform Wo f a has an abrupt transition the minimization of can yield a smoother solution W h a which oscillat

Download Pdf Manuals

image

Related Search

Related Contents

Manual PDF  TM8100 User`s Guide    Manual de Instalacion del Aplicativo OANET  Aclarando conceptos: mcm y MCD - Recursos  取扱説明書  ATAVRMC200 Hardware User Guide  DCR-DVD450E  I32_UserGuide_English  manuel - Cabasse  

Copyright © All rights reserved.
Failed to retrieve file