提升小波变换
简单实现了一下提升小波变换。包括Haar，（1，1），5/3，9/7 提升结构。欢迎有兴趣的同学一起讨论
function =lifting_transform_forward(im,levels,wavetype)
wavelet = wavetype;
for k =1:levels
=Nonlinear_lifting_forward(im,wavelet);
= Nonlinear_lifting_forward(L',wavelet);
= Nonlinear_lifting_forward(H',wavelet);
im = LL';
end
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 = lifting_im';
function =lifting_transform_inverse(im,levels,wavetype)
wavelet = wavetype;
= size(im);
min_len = m/2^levels;
LL = im(,);
for k= levels:-1:1
im_len = m/2^k;
LH{k} = im(,);
HL{k} = im(,);
HH{k} = im(,);
end
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';
function =Nonlinear_lifting_forward(im,wavelet)
= size(im);
%%%%split
for i=1:2:(width-1)
Xo((i+1)/2,:) = im(i,:);
end
for i=2:2:width
Xe(i/2,:) = im(i,:);
end
=size(Xe);
%%%% 5/3 lifting %%%%%%%%%%%%%%
if isequal(wavelet, '5.3')
next_Xe(,:) = Xe(,:);
next_Xe(sub_im_height,:) = Xe(sub_im_height,:);
d(:,:) = Xo(:,:) - (0.5*Xe(:,:)+0.5*next_Xe(:,:));
up_d(,:) = d(,:);
up_d(1,:) = d(1,:);
c(:,:) = Xe(:,:) +(0.25*(up_d(:,:)+d(:,:)));
end
%%%% 9/7 lifting %%%%%%%%%%%%%%
if isequal(wavelet, '9.7')
alfa = -1.586134342;
beta = -0.052980118;
gamma =0.882911075;
delta =0.443506852;
k =1.149604398;
C0 = Xe;
D0 = Xo;
C0_1(,:) = C0(,:);
C0_1(sub_im_height,:) = C0(sub_im_height,:);
D1(:,:) = D0(:,:)+alfa.*(C0(:,:)+C0_1(:,:));
D1_1(,:) = D1(,:);
D1_1(1,:) = D1(1,:);
C1(:,:) = C0(:,:)+beta.*(D1(:,:)+D1_1(:,:));
C1_1(,:) = C1(,:);
C1_1(sub_im_height,:) = C1(sub_im_height,:);
D2(:,:) = D1(:,:)+gamma.*(C1(:,:)+C1_1(:,:));
D2_1(,:) = D2(,:);
D2_1(1,:) = D2(1,:);
C2(:,:) = C1(:,:)+delta.*(D2(:,:)+D2_1(:,:));
c = k.* C2;
d = D2./k;
end
ifisequal(wavelet, '1.1')
c(:,:) = 0.5*(Xe(:,:)+Xo(:,:));
d(:,:) = Xo(:,:)-c(:,:);
end
if isequal(wavelet,'Haar')
c0 = Xe;
d0 = Xo;
d(:,:) = d0(:,:)-c0(:,:);
c(:,:) = c0(:,:)+0.5.*d(:,:);
end
% figure;
% imshow(uint8(c));
% figure;
% imshow(uint8(d));
%
function =Nonlinear_lifting_inverse(L,H,wavelet)
c = L;
d = H;
=size(L);
%%%% 5/3 lifting %%%%%%%%%%%%%%
if isequal(wavelet, '5.3')
up_d(,:) = d(,:);
up_d(1,:) = d(1,:);
Xe(:,:) = c(:,:) - (0.25*(up_d(:,:)+d(:,:)));
next_Xe(,:) = Xe(,:);
next_Xe(sub_im_height,:) = Xe(sub_im_height,:);
Xo(:,:) = d(:,:) + (0.5*Xe(:,:)+0.5*next_Xe(:,:));
end
%%%% 9/7 lifting %%%%%%%%%%%%%%
if isequal(wavelet, '9.7')
alfa = -1.586134342;
beta = -0.052980118;
gamma =0.882911075;
delta =0.443506852;
k =1.149604398;
D2 = d.*k;
C2 = c./k;
D2_1(,:) = D2(,:);
D2_1(1,:) = D2(1,:);
C1(:,:) = C2(:,:)-delta.*(D2(:,:)+D2_1(:,:));
C1_1(,:) = C1(,:);
C1_1(sub_im_height,:) = C1(sub_im_height,:);
D1(:,:) = D2(:,:)-gamma.*(C1(:,:)+C1_1(:,:));
D1_1(,:) = D1(,:);
D1_1(1,:) = D1(1,:);
C0(:,:) = C1(:,:)-beta.*(D1(:,:)+D1_1(:,:));
C0_1(,:) = C0(,:);
C0_1(sub_im_height,:) = C0(sub_im_height,:);
D0(:,:) = D1(:,:)-alfa.*(C0(:,:)+C0_1(:,:));
Xe = C0;
Xo = D0;
end
ifisequal(wavelet, '1.1')
d0(:,:) = d(:,:)+c(:,:);
c0(:,:) = 2.*(c(:,:)-0.5.*d0(:,:));
Xe = c0;
Xo = d0;
end
if isequal(wavelet,'Haar')
c0(:,:) = c(:,:)-0.5.*d(:,:);
d0(:,:) = d(:,:)+c0(:,:);
Xe = c0;
Xo = d0;
end
res(,:) = Xo(:,:);
res(,:) = Xe(:,:);
function testing_lifting_linjie
close all;
clear all;
clc;
im = double(imread('lena.bmp'));
figure;
subplot(221);
imshow(uint8(im));
title('Original Image');
levels =3;
= size(im);
lifting_im = lifting_transform_forward(im,levels,'5.3');
subplot(222);
imshow(uint8(lifting_im));
title('lifting coefficients');
result_im = lifting_transform_inverse(lifting_im,levels,'5.3');
high =sum(sum(abs(lifting_im(:,:))))
MSE=sum(sum((result_im-im).^2))/(im_h*im_w);
PSNR = 20*log((255*255)/MSE);
disp([' PSNR =',num2str(PSNR)]);
subplot(223);
imshow(uint8(result_im));
title('reconstruction image') 顶啊!!
很好的东西！
很好的东西！希望能学习学习！
西安电子科技大学的同胞！
我是西安的，可以联系吗？ 谢谢 学到东西了 谢谢，赞～ 好东西，谢谢版主，辛苦了
谢
:handshake ,:victory: 我下载了代码.在使用过程中我发现了一个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 编辑 ]
回复 8楼 zhu1982lin 的帖子
我按照你的改了 可还是不行啊！一直提示有错！??? Error using ==> cat
CAT arguments dimensions are not consistent.
Error in ==> lifting_transform_forward at 13
lifting_im = cat(1,temp1',temp2'); 亲有木有提升小波包的代码哇 迷失的傻 发表于 2012-9-12 15:18 static/image/common/back.gif
亲有木有提升小波包的代码哇
function varargout = lwt(x,LS,varargin)
%LWT Lifting wavelet decomposition 1-D.
% LWT performs a 1-D lifting wavelet decomposition
% with respect to a particular lifted wavelet that you specify.
%
% = LWT(X,W) computes the approximation
% coefficients vector CA and detail coefficients vector CD,
% obtained by a lifting wavelet decomposition, of
% the vector X. W is a lifted wavelet name (see LIFTWAVE).
%
% X_InPlace = LWT(X,W) computes the approximation and
% detail coefficients. These coefficients are stored in-place:
% CA = X_InPlace(1:2:end) and CD = X_InPlace(2:2:end)
%
% LWT(X,W,LEVEL) computes the lifting wavelet decomposition
% at level LEVEL.
%
% X_InPlace = LWT(X,W,LEVEL,'typeDEC',typeDEC) or
% = LWT(X,W,LEVEL,'typeDEC',typeDEC) with
% typeDEC = 'w' or 'wp' computes the wavelet or the
% wavelet packet decomposition using lifting, at level LEVEL.
%
% Instead of a lifted wavelet name, you may use the associated
% lifting scheme LS:
% LWT(X,LS,...) instead of LWT(X,W,...).
%
% For more information about lifting schemes type: lsinfo.
%
% See also ILWT.
% M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 02-Feb-2000.
% Last Revision: 10-Jul-2003.
% Copyright 1995-2004 The MathWorks, Inc.
% $Revision: 1.1.6.3 $$Date: 2004/04/13 00:39:57 $
% Check arguments.
nbIn = nargin;
switch nbIn
case {2,3,5} ,
case {0,1} , error('Not enough input arguments.');
case {4} , error('Invalid number of input arguments.');
otherwise, error('Too many input arguments.');
end
if nargout>3
error('Too many output arguments.');
end
% Default: level and typeDEC.
level = 1;
typeDEC = 'w';
if nargin>2
level = varargin{1};
for k = 2:2:length(varargin)-1
argName = lower( varargin{k});
switch argName
case 'typedec' , typeDEC = varargin{k+1};
end
end
end
% if level>1
% msg = nargoutchk(0,1,nargout); error(msg);
% end
if ischar(LS) , LS = liftwave(LS); end
%===================%
% LIFTING ALGORITHM %
%===================%
% Splitting.
lx = length(x);
firstIdxAPP = 1;
firstIdxDET = 1+mod(firstIdxAPP,2);
idxAPP = firstIdxAPP:2:lx;
idxDET = firstIdxDET:2:lx;
lenAPP = length(idxAPP);
lenDET = length(idxDET);
% Lifting.
NBL = size(LS,1);
LStype = LS{NBL,3};
for k = 1:NBL-1
liftTYPE = LS{k,1};
liftFILT = LS{k,2};
DF = LS{k,3};
switch liftTYPE
case 'p' ,
x(idxAPP) = x(idxAPP) + ...
lsupdate('v',x(idxDET),liftFILT,DF,lenAPP,LStype);
case 'd' ,
x(idxDET) = x(idxDET) + ...
lsupdate('v',x(idxAPP),liftFILT,DF,lenDET,LStype);
end
end
% Normalization.
if isempty(LStype)
x(idxAPP) = LS{NBL,1}*x(idxAPP);
x(idxDET) = LS{NBL,2}*x(idxDET);
end
%========================================================================%
% Recursion if level > 1.
if level>1
x(idxAPP) = lwt(x(idxAPP),LS,level-1,'typeDEC',typeDEC);
if isequal(typeDEC,'wp')
x(idxDET) = lwt(x(idxDET),LS,level-1,'typeDEC',typeDEC);
end
end
switch nargout
case 1 , varargout = {x};
case 2 , varargout = {x(idxAPP),x(idxDET)};
case 3 , varargout = {x,x(idxAPP),x(idxDET)};
end
页:
[1]