|
楼主 |
发表于 2006-3-31 08:36
|
显示全部楼层
提升法97经典程序
- %% 本程序实现任意偶数大小图像第二代双正交97提升小波变换 %% 注1: 采用标准正交方法,对行列采用不同矩阵(和matlab里不同)
- %% 注2: 为了保证正交,所有边界处理,全部采用循环处理
- %% 注3: 正交性验证,将单位阵带入函数,输出仍是单位阵(matlab不具有此性质)
- %% 注4: 此程序是矩阵实现,所以图像水平分量和垂直分量估计被交换位置
- %% 注5: 此程序实现的是类小波(wavelet-like)变换,是介于小波包变换与小波变换之间的变换
- %% 注6: 此程序每层变换相对原图像矩阵,产生的矩阵都是正交阵,这和小波包一致
- %% 注7: 但小波变换每层产生的矩阵,是相对每个待分解子块的正交矩阵,而不是原图像的正交矩阵
- %% 注8: 且小波变换产生的正交矩阵维数,随分解层数2分减少
- %% 注9: 提升系数可以在MATLAB7.0以上版本,用liftwave('9.7')获取,这里直接给出,考虑兼容性
- %% 注10:由于MATLAB数组下标从1开始,所以注意奇偶序列的变化
- %% 注11:d为对偶上升,即预测;p为原上升,即更新 %% 编程人 沙威 安徽大学
- %% 编程时间 2004年12月18日 %% x输入图像,y输出图像
- %% flag_trans为正变换或反变换标志,0执行正变换,1执行反变换
- %% flag_max,是否最大层数变换标志,0执行用户设定层数,1执行最大层数变换
- %% layer,用户层数设置(小于最大层) function y=db97(x,flag_trans,flag_max,layer); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % 1.输入参数检查 % 矩阵维数判断
- [sa,sb]=size(x); if (sa~=sb) % 防止非图像数据
- errordlg('非图像数据!');
- error('非图像数据!');
- end; % 变换标志判断
- [sa,sb]=size(flag_trans);
- if ((sa~=1) | (sb~=1)) % 变换标志错误
- errordlg('变换标志错误!');
- error('变换标志错误!');
- end; if ((flag_trans~=1) & (flag_trans~=0)) % 变换标志错误
- errordlg('变换标志错误!');
- error('变换标志错误!');
- end; % 最大层数标志判断
- [sa,sb]=size(flag_max);
- if ((sa~=1) | (sb~=1)) % 最大层数标志错误
- errordlg('最大层数标志错误!');
- error('最大层数标志错误!');
- end; if ((flag_max~=1) & (flag_max~=0)) % 最大层数标志错误
- errordlg('最大层数标志错误!');
- error('最大层数标志错误!');
- end; % 用户设置层数判断
- if (flag_max~=1) [sa,sb]=size(layer);
- if ((sa~=1) | (sb~=1)) % 层数设置错误
- errordlg('层数设置错误!');
- error('层数设置错误!');
- end; if (flag_max<0) % 层数设置错误
- errordlg('层数设置错误!');
- error('层数设置错误!');
- end;
- end;
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % 2.提升系数确定
- % t1=liftwave('9.7'); % 获取提升系数(MATLAB7.0以后) d1=[-1.586100000000000e+000,-1.586134342069360e+000];
- p1=[1.079600000000000e+000,-5.298011857188560e-002];
- d2=[-8.829110755411875e-001,-8.829110755411875e-001];
- p2=[4.435068520511142e-001,1.576123746148364e+000];
- d3=-8.698644516247808e-001;
- p3=-1.149604398860242e+000;
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % 3.分解层数确定
- % 采用用户输入和自动给出最大层数两种方法 N=length(x); % 矩阵大小
- S=N; % 变量
- s=log2(N); % 最大循环次数
- n1=N/2; % 初始一半矩阵大小
- n2=N; % 初始矩阵大小
- u=0; % 初始值 % 对非2的整数幂大小图像确定最大分解层数
- for ss=1:s
- if (mod(S,2)==0)
- u=u+1;
- S=S/2;
- end;
- end;
- u=u-1; % 分解最大层数减1(后面的边界处理造成) % 最大层数确定
- if (flag_max==0) % 手动输入
- T=layer; % 用户输入值
- else % 自动确定最大层数
- T=u; % 分解最大层数
- end
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % 4.最大层数和图像大小检查 if (T>u) % 防止用户层数越界
- errordlg('已超过最大分解层数!或者非偶数大小图像!');
- error('已超过最大分解层数!或者非偶数大小图像!');
- end; if (mod(N,2)~=0) % 防止图像大小错误
- errordlg('非偶数大小图像!');
- error('非偶数大小图像!');
- end;
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % 5.提升法正变换 if (flag_trans==0)
- for time=1:T; % 行正变换
- % d;
- x1(n1,:)=x(n2,:)+d1(2)*x(n2-1,:)+d1(1)*x(1,:);
- x1([1:n1-1],:)=x([2:2:n2-2],:)+d1(2)*x([1:2:n2-3],:)+d1(1)*x([3:2:n2-1],:);
- % p;
- x(1,:)=x(1,:)+p1(2)*x1(n1,:)+p1(1)*x1(1,:);
- x([2:n1],:)=x([3:2:n2-1],:)+p1(2)*x1([1:n1-1],:)+p1(1)*x1([2:n1],:);
- x([n1+1:n2],:)=x1([1:n1],:);
- % d;
- x(n1+1,:)=x(n1+1,:)+d2(2)*x(n1,:)+d2(1)*x(1,:);
- x([n1+2:n2],:)=x([n1+2:n2],:)+d2(2)*x([1:n1-1],:)+d2(1)*x([2:n1],:);
- % p;
- x(n1,:)=x(n1,:)+p2(2)*x(n1+1,:)+p2(1)*x(n1+2,:);
- x(n1-1,:)=x(n1-1,:)+p2(2)*x(n2,:)+p2(1)*x(n1+1,:);
- x([1:n1-2],:)=x([1:n1-2],:)+p2(2)*x([n1+2:n2-1],:)+p2(1)*x([n1+3:n2],:);
- % 归一
- x([1:n1],:)=p3*x([1:n1],:);
- x([n1+1:n2],:)=d3*x([n1+1:n2],:); clear x1;
- % 列正变换
- % d;
- x1(:,[1:n1])=x(:,[2:2:n2]);
- % p;
- x(:,1)=x(:,1)-d1(1)*x1(:,n1)-d1(2)*x1(:,1);
- x(:,[2:n1])=x(:,[3:2:n2-1])-d1(1)*x1(:,[1:n1-1])-d1(2)*x1(:,[2:n1]);
- x(:,[n1+1:n2])=x1(:,[1:n1]);
- % d;
- x(:,n2)=x(:,n2)-p1(1)*x(:,n1)-p1(2)*x(:,1);
- x(:,[n1+1:n2-1])=x(:,[n1+1:n2-1])-p1(1)*x(:,[1:n1-1])-p1(2)*x(:,[2:n1]);
- % p;
- x(:,n1,:)=x(:,n1)-d2(1)*x(:,n2)-d2(2)*x(:,n1+1);
- x(:,[1:n1-1])=x(:,[1:n1-1])-d2(1)*x(:,[n1+1:n2-1])-d2(2)*x(:,[n1+2:n2]);
- % d;
- x(:,n1+1)=x(:,n1+1)-p2(1)*x(:,n1-1)-p2(2)*x(:,n1);
- x(:,n1+2)=x(:,n1+2)-p2(1)*x(:,n1)-p2(2)*x(:,1);
- x(:,[n1+3:n2])=x(:,[n1+3:n2])-p2(1)*x(:,[1:n1-2])-p2(2)*x(:,[2:n1-1]);
- % 归一
- x(:,[1:n1])=d3*x(:,[1:n1]);
- x(:,[n1+1:n2])=p3*x(:,[n1+1:n2]); clear x1;
- n2=n2/2; % 原大小
- n1=n2/2; % 一半大小
- end;
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % 6.提升法反变换 else
- n2=N/(2.^(T-1)); % 分解最小子块维数
- n1=n2/2;
- for time=1:T; % 行反变换
- % 去归一
- x([1:n1],:)=x([1:n1],:)/p3;
- x([n1+1:n2],:)=x([n1+1:n2],:)/d3; % 反p;
- x(n1,:)=x(n1,:)-p2(2)*x(n1+1,:)-p2(1)*x(n1+2,:);
- x(n1-1,:)=x(n1-1,:)-p2(2)*x(n2,:)-p2(1)*x(n1+1,:);
- x([1:n1-2],:)=x([1:n1-2],:)-p2(2)*x([n1+2:n2-1],:)-p2(1)*x([n1+3:n2],:);
- % 反d;
- x(n1+1,:)=x(n1+1,:)-d2(2)*x(n1,:)-d2(1)*x(1,:);
- x([n1+2:n2],:)=x([n1+2:n2],:)-d2(2)*x([1:n1-1],:)-d2(1)*x([2:n1],:);
- % 反p;
- x1(1,:)=x(1,:)-p1(2)*x(n2,:)-p1(1)*x(n1+1,:);
- x1([2:n1],:)=x([2:n1],:)-p1(2)*x([n1+1:n2-1],:)-p1(1)*x([n1+2:n2],:);
- % 反d;
- x(n2,:)=x(n2,:)-d1(2)*x1(n1,:)-d1(1)*x1(1,:);
- x([2:2:n2-2],:)=x([n1+1:n2-1],:)-d1(2)*x1([1:n1-1],:)-d1(1)*x1([2:n1],:);
- % 偶数
- x([1:2:n2-1],:)=x1([1:n1],:);
- clear x1;
- % 列反变换
- % 归一
- x(:,[1:n1])=x(:,[1:n1])/d3;
- x(:,[n1+1:n2])=x(:,[n1+1:n2])/p3; % 反d;
- x(:,n1+1)=x(:,n1+1)+p2(1)*x(:,n1-1)+p2(2)*x(:,n1);
- x(:,n1+2)=x(:,n1+2)+p2(1)*x(:,n1)+p2(2)*x(:,1);
- x(:,[n1+3:n2])=x(:,[n1+3:n2])+p2(1)*x(:,[1:n1-2])+p2(2)*x(:,[2:n1-1]);
- % 反p;
- x(:,n1,:)=x(:,n1)+d2(1)*x(:,n2)+d2(2)*x(:,n1+1);
- x(:,[1:n1-1])=x(:,[1:n1-1])+d2(1)*x(:,[n1+1:n2-1])+d2(2)*x(:,[n1+2:n2]);
- % 反d;
- x(:,n2)=x(:,n2)+p1(1)*x(:,n1)+p1(2)*x(:,1);
- x(:,[n1+1:n2-1])=x(:,[n1+1:n2-1])+p1(1)*x(:,[1:n1-1])+p1(2)*x(:,[2:n1]);
- % 反p;
- x1(:,1)=x(:,1)+d1(1)*x(:,n2)+d1(2)*x(:,n1+1);
- x1(:,[2:n1])=x(:,[2:n1])+d1(1)*x(:,[n1+1:n2-1])+d1(2)*x(:,[n1+2:n2]); % 奇偶
- x(:,[2:2:n2])=x(:,[n1+1:n2]);
- x(:,[1:2:n2-1])=x1(:,[1:n1]); clear x1;
- n2=n2*2; % 原大小
- n1=n2/2; % 一半大小 end;
- end;
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % 7.结果输出 y=x;
- % 传输最后结果 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % 8.内存清理 clear x;
- clear flag_max;
- clear layer;
- clear flag_trans;
- clear N;
- clear n1;
- clear n2;
- clear s;
- clear ss;
- clear u;
- clear d1;
- clear d2;
- clear d3;
- clear p1;
- clear p2;
- clear p3;
- clear sa;
- clear sb;
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
复制代码
[ 本帖最后由 yejet 于 2006-8-31 20:32 编辑 ] |
|