提升法97经典程序.txt

本程序实现任意偶数大小图像第二代双正交97提升小波变换

提升法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为原上升，即更新
%%  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;