提升小波变换

lanpijiu

简单实现了一下提升小波变换。包括Haar，（1，1），5/3，9/7 提升结构。
欢迎有兴趣的同学一起讨论

1. function [res]=lifting_transform_forward(im,levels,wavetype)
   
2. wavelet = wavetype;
   
3. for k =1:levels
   
4.     [L,H]=Nonlinear_lifting_forward(im,wavelet);
   
5.     [LL,LH{k}] = Nonlinear_lifting_forward(L',wavelet);
   
6.     [HL{k},HH{k}] = Nonlinear_lifting_forward(H',wavelet);
   
7.     im = LL';
   
8. end
   
9. 
   
10. for k= levels:-1:1
    
11.     temp1 =cat(2,LL,LH{k});
    
12.     temp2 =cat(2,HL{k},HH{k});
    
13.     lifting_im = cat(1,temp1,temp2);
    
14.     LL = lifting_im;
    
15. end
    
16. res = lifting_im';

复制代码

1. function [res]=lifting_transform_inverse(im,levels,wavetype)
   
2. wavelet = wavetype;
   
3. [m,n] = size(im);
   
4. min_len = m/2^levels;
   
5. LL = im([1:min_len],[1:min_len]);
   
6. for k= levels:-1:1
   
7.     im_len = m/2^k;
   
8.     LH{k} = im([1:im_len],[im_len+1:im_len*2]);
   
9.     HL{k} = im([im_len+1:im_len*2],[1:im_len]);
   
10.     HH{k} = im([im_len+1:im_len*2],[im_len+1:im_len*2]);
    
11. end
    
12. for k =levels:-1:1
    
13.     up = Nonlinear_lifting_inverse(LL,LH{k},wavelet);
    
14.     down = Nonlinear_lifting_inverse(HL{k},HH{k},wavelet);
    
15.     im  = Nonlinear_lifting_inverse(up',down',wavelet);
    
16.     LL= im';
    
17. end
    
18. res = im';

复制代码

1. function [c,d]=Nonlinear_lifting_forward(im,wavelet)
   
2. [width,height] = size(im);
   
3. %%%%split
   
4. for i=1:2:(width-1)
   
5.     Xo((i+1)/2,:) = im(i,:);
   
6. end
   
7. for i=2:2:width
   
8.     Xe(i/2,:) = im(i,:);
   
9. end
   
10. [sub_im_height,sub_im_width]=size(Xe);
    
11. %%%%   5/3 lifting %%%%%%%%%%%%%%
    
12. if isequal(wavelet, '5.3')
    
13.     next_Xe([1:sub_im_height-1],:) = Xe([2:sub_im_height],:);
    
14.     next_Xe(sub_im_height,:) = Xe(sub_im_height,:);
    
15.     d(:,:) = Xo(:,:) - (0.5*Xe(:,:)+0.5*next_Xe(:,:));
    
16.     up_d([2:sub_im_height],:) = d([1:sub_im_height-1],:);
    
17.     up_d(1,:) = d(1,:);
    
18.     c(:,:) = Xe(:,:) +(0.25*(up_d(:,:)+d(:,:)));
    
19. end
    
20. %%%%   9/7 lifting %%%%%%%%%%%%%%
    
21. if isequal(wavelet, '9.7')
    
22.     alfa    = -1.586134342;
    
23.     beta    = -0.052980118;
    
24.     gamma   =  0.882911075;
    
25.     delta   =  0.443506852;
    
26.     k       =  1.149604398;
    
27.     C0 = Xe;
    
28.     D0 = Xo;
    
29.     C0_1([1:sub_im_height-1],:) = C0([2:sub_im_height],:);
    
30.     C0_1(sub_im_height,:) = C0(sub_im_height,:);
    
31.     D1(:,:) = D0(:,:)+alfa.*(C0(:,:)+C0_1(:,:));
    
32.     D1_1([2:sub_im_height],:) = D1([1:sub_im_height-1],:);
    
33.     D1_1(1,:) = D1(1,:);
    
34.     C1(:,:) = C0(:,:)+beta.*(D1(:,:)+D1_1(:,:));
    
35.     C1_1([1:sub_im_height-1],:) = C1([2:sub_im_height],:);
    
36.     C1_1(sub_im_height,:) = C1(sub_im_height,:);
    
37.     D2(:,:) = D1(:,:)+gamma.*(C1(:,:)+C1_1(:,:));
    
38.     D2_1([2:sub_im_height],:) = D2([1:sub_im_height-1],:);
    
39.     D2_1(1,:) = D2(1,:);
    
40.     C2(:,:) = C1(:,:)+delta.*(D2(:,:)+D2_1(:,:));
    
41.     c = k.* C2;
    
42.     d = D2./k;
    
43. end
    
44. if  isequal(wavelet, '1.1')
    
45.     c(:,:) = 0.5*(Xe(:,:)+Xo(:,:));
    
46.     d(:,:) = Xo(:,:)-c(:,:);
    
47. end
    
48. 
    
49. if isequal(wavelet,'Haar')
    
50.     c0 = Xe;
    
51.     d0 = Xo;
    
52.     d(:,:) = d0(:,:)-c0(:,:);
    
53.     c(:,:) = c0(:,:)+0.5.*d(:,:);
    
54. end
    
55. % figure;
    
56. % imshow(uint8(c));
    
57. % figure;
    
58. % imshow(uint8(d));
    
59. %

复制代码

1. function [res]=Nonlinear_lifting_inverse(L,H,wavelet)
   
2. c = L;
   
3. d = H;
   
4. [sub_im_height,sub_im_width]=size(L);
   
5. %%%%   5/3 lifting %%%%%%%%%%%%%%
   
6. if isequal(wavelet, '5.3')
   
7.     up_d([2:sub_im_height],:) = d([1:sub_im_height-1],:);
   
8.     up_d(1,:) = d(1,:);
   
9.     Xe(:,:) = c(:,:) - (0.25*(up_d(:,:)+d(:,:)));
   
10.    
    
11.     next_Xe([1:sub_im_height-1],:) = Xe([2:sub_im_height],:);
    
12.     next_Xe(sub_im_height,:) = Xe(sub_im_height,:);
    
13.     Xo(:,:) = d(:,:) + (0.5*Xe(:,:)+0.5*next_Xe(:,:));
    
14. end
    
15. %%%%   9/7 lifting %%%%%%%%%%%%%%
    
16. if isequal(wavelet, '9.7')
    
17.     alfa    = -1.586134342;
    
18.     beta    = -0.052980118;
    
19.     gamma   =  0.882911075;
    
20.     delta   =  0.443506852;
    
21.     k       =  1.149604398;
    
22.     D2 = d.*k;
    
23.     C2 = c./k;
    
24.     D2_1([2:sub_im_height],:) = D2([1:sub_im_height-1],:);
    
25.     D2_1(1,:) = D2(1,:);
    
26.     C1(:,:) = C2(:,:)-delta.*(D2(:,:)+D2_1(:,:));
    
27.     C1_1([1:sub_im_height-1],:) = C1([2:sub_im_height],:);
    
28.     C1_1(sub_im_height,:) = C1(sub_im_height,:);
    
29.     D1(:,:) = D2(:,:)-gamma.*(C1(:,:)+C1_1(:,:));
    
30.     D1_1([2:sub_im_height],:) = D1([1:sub_im_height-1],:);
    
31.     D1_1(1,:) = D1(1,:);
    
32.     C0(:,:) = C1(:,:)-beta.*(D1(:,:)+D1_1(:,:));
    
33.     C0_1([1:sub_im_height-1],:) = C0([2:sub_im_height],:);
    
34.     C0_1(sub_im_height,:) = C0(sub_im_height,:);
    
35.     D0(:,:) = D1(:,:)-alfa.*(C0(:,:)+C0_1(:,:));
    
36.      Xe = C0;
    
37.      Xo = D0;
    
38. end
    
39. if  isequal(wavelet, '1.1')
    
40.     d0(:,:) = d(:,:)+c(:,:);
    
41.     c0(:,:) = 2.*(c(:,:)-0.5.*d0(:,:));
    
42.      Xe = c0;
    
43.      Xo = d0;
    
44. end
    
45. if isequal(wavelet,'Haar')
    
46.     c0(:,:) = c(:,:)-0.5.*d(:,:);
    
47.     d0(:,:) = d(:,:)+c0(:,:);
    
48.     Xe = c0;
    
49.     Xo = d0;
    
50. end
    
51. res([1:2:sub_im_height*2-1],:) = Xo(:,:);
    
52. res([2:2:sub_im_height*2],:) = Xe(:,:);

复制代码

1. function testing_lifting_linjie
   
2. close all;
   
3. clear all;
   
4. clc;
   
5. 
   
6. im = double(imread('lena.bmp'));
   
7. figure;
   
8. subplot(221);
   
9. imshow(uint8(im));
   
10. title('Original Image');  
    
11. 
    
12. levels =3;
    
13. [im_h,im_w] = size(im);
    
14. lifting_im = lifting_transform_forward(im,levels,'5.3');
    
15. subplot(222);
    
16. imshow(uint8(lifting_im));
    
17. title('lifting coefficients');  
    
18. 
    
19. result_im = lifting_transform_inverse(lifting_im,levels,'5.3');
    
20. high =  sum(sum(abs(lifting_im(:,:))))
    
21. MSE=sum(sum((result_im-im).^2))/(im_h*im_w);
    
22. PSNR = 20*log((255*255)/MSE);
    
23. disp(['     PSNR =  ',num2str(PSNR)]);
    
24. subplot(223);
    
25. imshow(uint8(result_im));
    
26. title('reconstruction image')

复制代码

pltlemc

顶啊!!

xingchen

很好的东西！
希望能学习学习！
西安电子科技大学的同胞！
我是西安的，可以联系吗？

loveming

谢谢 学到东西了

davlee

谢谢，赞～

zongkui

好东西，谢谢版主，辛苦了

祥云

:handshake ,:victory:

zhu1982lin

我下载了代码.
在使用过程中我发现了一个bug.
由于我只需要5.3提升小波,其他的就没去检查了.
在5.3提升小波的实现中,以下代码应该存在bug:

for k= levels:-1:1
    temp1 =cat(2,LL,LH{k});
    temp2 =cat(2,HL{k},HH{k});
    lifting_im = cat(1,temp1,temp2);
    LL = lifting_im;
end
res = im';

和
for k =levels:-1:1
    up = Nonlinear_lifting_inverse(LL,LH{k},wavelet);
    down = Nonlinear_lifting_inverse(HL{k},HH{k},wavelet);
    im  = Nonlinear_lifting_inverse(up',down',wavelet);
    LL= im';
end

应该改为:
for k= levels:-1:1
    temp1 =cat(1,LL,LH{k});
    temp2 =cat(1,HL{k},HH{k});
    lifting_im = cat(1, temp1' , temp2' );
    LL = lifting_im';
end
res = lifting_im;

和
for k =levels:-1:1
    up = Nonlinear_lifting_inverse(LL' , LH{k}' , wavelet);
    down = Nonlinear_lifting_inverse(HL{k}' , HH{k}' , wavelet);
    im  = Nonlinear_lifting_inverse(up', down', wavelet);
    LL= im;
end
res = im;

我还修改了函数 Nonlinear_lifting_forward 和 Nonlinear_lifting_inverse 中 5.3提升小波的计算方法,主要是循环扩展,和 Y(2n) = X(2n) + [(Y(2n-1) + Y(2n+1) + 2)/4]; (原文使用的公式为:Y(2n) = X(2n) + [(Y(2n-1) + Y(2n+1))/4])

[ 本帖最后由 zhu1982lin 于 2009-1-29 05:30 编辑 ]

qqvirile

我按照你的改了 可还是不行啊！一直提示有错！
??? Error using ==> cat
CAT arguments dimensions are not consistent.

Error in ==> lifting_transform_forward at 13
    lifting_im = cat(1,temp1',temp2');

迷失的傻

亲有木有提升小波包的代码哇

funi

迷失的傻 发表于 2012-9-12 15:18

                               
登录/注册后可看大图

亲有木有提升小波包的代码哇

1. 
   
2. function varargout = lwt(x,LS,varargin)
   
3. %LWT Lifting wavelet decomposition 1-D.
   
4. %   LWT performs a 1-D lifting wavelet decomposition
   
5. %   with respect to a particular lifted wavelet that you specify.
   
6. %
   
7. %   [CA,CD] = LWT(X,W) computes the approximation
   
8. %   coefficients vector CA and detail coefficients vector CD,
   
9. %   obtained by a lifting wavelet decomposition, of
   
10. %   the vector X. W is a lifted wavelet name (see LIFTWAVE).
    
11. %
    
12. %   X_InPlace = LWT(X,W) computes the approximation and
    
13. %   detail coefficients. These coefficients are stored in-place:
    
14. %     CA = X_InPlace(1:2:end) and CD = X_InPlace(2:2:end)
    
15. %
    
16. %   LWT(X,W,LEVEL) computes the lifting wavelet decomposition
    
17. %   at level LEVEL.
    
18. %
    
19. %   X_InPlace = LWT(X,W,LEVEL,'typeDEC',typeDEC) or
    
20. %   [CA,CD] = LWT(X,W,LEVEL,'typeDEC',typeDEC) with
    
21. %   typeDEC = 'w' or 'wp' computes the wavelet or the
    
22. %   wavelet packet decomposition using lifting, at level LEVEL.
    
23. %
    
24. %   Instead of a lifted wavelet name, you may use the associated
    
25. %   lifting scheme LS:
    
26. %     LWT(X,LS,...) instead of LWT(X,W,...).
    
27. %
    
28. %   For more information about lifting schemes type: lsinfo.
    
29. %
    
30. %   See also ILWT.
    
31. 
    
32. %   M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 02-Feb-2000.
    
33. %   Last Revision: 10-Jul-2003.
    
34. %   Copyright 1995-2004 The MathWorks, Inc.
    
35. %   $Revision: 1.1.6.3 $  $Date: 2004/04/13 00:39:57 $
    
36. 
    
37. % Check arguments.
    
38. nbIn = nargin;
    
39. switch nbIn
    
40.     case {2,3,5} ,
    
41.     case {0,1} , error('Not enough input arguments.');
    
42.     case {4}   , error('Invalid number of input arguments.');
    
43.     otherwise  , error('Too many input arguments.');
    
44. end
    
45. if nargout>3
    
46.     error('Too many output arguments.');
    
47. end
    
48. 
    
49. % Default: level and typeDEC.
    
50. level = 1;
    
51. typeDEC = 'w';
    
52. if nargin>2
    
53.     level = varargin{1};
    
54.     for k = 2:2:length(varargin)-1
    
55.         argName = lower( varargin{k});
    
56.         switch argName
    
57.             case 'typedec' , typeDEC = varargin{k+1};
    
58.         end
    
59.     end
    
60. end
    
61. % if level>1
    
62. %     msg = nargoutchk(0,1,nargout); error(msg);
    
63. % end
    
64. if ischar(LS) , LS = liftwave(LS); end
    
65. 
    
66. %===================%
    
67. % LIFTING ALGORITHM %
    
68. %===================%
    
69. % Splitting.
    
70. lx = length(x);
    
71. firstIdxAPP = 1;
    
72. firstIdxDET = 1+mod(firstIdxAPP,2);
    
73. idxAPP = firstIdxAPP:2:lx;
    
74. idxDET = firstIdxDET:2:lx;
    
75. lenAPP = length(idxAPP);
    
76. lenDET = length(idxDET);
    
77. 
    
78. % Lifting.
    
79. NBL = size(LS,1);
    
80. LStype = LS{NBL,3};
    
81. for k = 1:NBL-1
    
82.     liftTYPE = LS{k,1};
    
83.     liftFILT = LS{k,2};
    
84.     DF       = LS{k,3};
    
85.     switch liftTYPE
    
86.        case 'p' ,
    
87.            x(idxAPP) = x(idxAPP) + ...
    
88.                lsupdate('v',x(idxDET),liftFILT,DF,lenAPP,LStype);
    
89.        case 'd' ,
    
90.            x(idxDET) = x(idxDET) + ...
    
91.                lsupdate('v',x(idxAPP),liftFILT,DF,lenDET,LStype);
    
92.     end
    
93. end
    
94. 
    
95. % Normalization.
    
96. if isempty(LStype)
    
97.     x(idxAPP) = LS{NBL,1}*x(idxAPP);
    
98.     x(idxDET) = LS{NBL,2}*x(idxDET);
    
99. end
    
100. %========================================================================%
     
101. 
     
102. % Recursion if level > 1.
     
103. if level>1
     
104.    x(idxAPP) = lwt(x(idxAPP),LS,level-1,'typeDEC',typeDEC);
     
105.    if isequal(typeDEC,'wp')
     
106.        x(idxDET) = lwt(x(idxDET),LS,level-1,'typeDEC',typeDEC);
     
107.    end
     
108. end
     
109. 
     
110. switch nargout
     
111.   case 1 , varargout = {x};
     
112.   case 2 , varargout = {x(idxAPP),x(idxDET)};
     
113.   case 3 , varargout = {x,x(idxAPP),x(idxDET)};
     
114. end
     

复制代码