1. 程式人生 > >影象的Mallat演算法分解(Matlab程式碼)

影象的Mallat演算法分解(Matlab程式碼)

Mallat 演算法的分析與綜合框架參考書上的資料很多,這裡就不多說了。

下面是我寫的關於影象的程式,分別是:一維分解,二維分解;一維合成,二維合成。最後是測試主程式。

謝謝參考,錯了請反饋一下!

%內部子函式,對一行(row)向量進行一次小波變換,利用fft實現
function y=mdec1(x,h,g)
%x 行陣列
%h為低通濾波器
%g為高通濾波器
%輸出: y 進行一級小波分解後的係數
lenx=size(x,2);                 %求行的長度
lenh=size(h,2);                 %求低通濾波器的長度
rh=h(end:-1:1);                 %rh是h的逆序  h0(n)=h(-n)
rrh=[zeros(1,(lenx-lenh)),rh];  %rrg在rg前插入(lenx-lenh))個零
rrh=circshift(rrh',1)';         %迴圈位移

rg=g(end:-1:1);                 %rg為高通濾波器的逆序g0(n)=g(-n)
rrg=[zeros(1,(lenx-lenh)),rg];  %rrg在rg前插入(lenx-lenh))個零
rrg=circshift(rrg',1)';         %迴圈位移

r1=dyaddown(ifft(fft(x).*fft(rrh,lenx)),1); %use para 1,下2奇取樣
r2=dyaddown(ifft(fft(x).*fft(rrg,lenx)),1); %use para 1,下2奇取樣
y=[r1,r2];                     %相加得到小波係數

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Y=mallatdec2(X,wname,level)
%輸入:X????? 載入的二維影象像數值;
%???? level? 小波分解次(級)數設定值(如果設定值超過最高可分解次數,按最高分解次數分解)
%????? wname? 小波名字wavelet name
%輸出:Y???? 多極小波分解後的小波係數矩陣
%wfilters:小波濾波器

[h,g]=wfilters(wname,'d');    %h,g分別為分解的低通和高通濾波器
X=double(X);
t=1;
hh=size(X,2);
Y=zeros(1,hh);
while t<=level   %迴圈
    %先進行行小波變換
    for row=1:hh
        Y(row,1:hh)=mdec1(X(row,1:hh),h,g) ;
    end
    %再進行列小波變換
    for col=1:hh
        temp=mdec1( Y(1:hh,col)',h,g); %Y(1:hh,col)為矩陣中的一列,轉置後為行向量,得到的結果temp也是行向量
        Y(1:hh,col)=temp';             %將行向量轉置為列向量儲存在Y中
    end
    t=t+1;                             %迴圈次數加1,大於level,迴圈終止
    hh=hh/2;                           %分解的矩陣的大小為原來的四分之一
    X=Y;                               %將Y值賦給X等待下一次分解                     
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%內部子函式,對一行小波係數進行重構
function y=mrec1(x,h,g)
%輸入:x 行陣列
%???? h為低通濾波器
%???? g為高通濾波器
%輸出: y 進行一級小波重構後值
lenx=size(x,2);

r3=dyadup(x(1,1:lenx*0.5),0);        %內插零use para 0,上取樣
r4=dyadup(x(1,(lenx*0.5+1):lenx),0); %內插零use para 0,上取樣
y=ifft(fft(r3,lenx).*fft(h,lenx))+ ifft(fft(r4,lenx).*fft(g,lenx));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%二維小波重構函式
function Y=mallatrec2(X,wname,level)
%輸入:X????? 載入的小波係數矩陣;
%?? ??level? 小波分解次(級)數設定值(如果設定值超過最高可分解次數,按最高分解次數分解)
%????? wname? 小波名字wavelet name
%輸出:Y???? 重構影象矩陣

[h,g]=wfilters(wname,'d');%h為重構的低通濾波器;g為重構高通濾波器

hz=size(X,2);
h1=hz/(2^(level-1));
while h1<=hz
    % 對列變換
    for col=1:h1
        temp=mrec1(X(1:h1,col)',h,g)';
        X(1:h1,col)=temp;
    end
    %再對行變換
    for row=1:h1
        temp=mrec1(X(row,1:h1),h,g);
        X(row,1:h1)=temp;
    end
    h1=h1*2;
end
Y=X;

 

%%%%%%%%%%%%%%%%%%%%%%%%%%

%測試函式(主函式)
clc;
clear;
close all;

X=imread('lena.bmp');;%路徑
figure(1);imshow(X);

X=double(X);
level=2;
A = mallatdec2(X,'db4',level);

figure(2);
imshow(uint8(A));title('多尺度分解影象');

Y= mallatrec2(A,'db4',level);

figure(3);
imshow(uint8(Y));title('重構影象');