%程序L7_2.m: 應(yīng)用小波變換與向量量化做圖像壓縮%
%利用LGB訓(xùn)練一層小波變換后四個(gè)頻帶的碼本%
function Wavelet_Vector_Comprssion()
%X=imread('j:\數(shù)字圖像處理－活用Matlab程序\CH7\lena.bmp','bmp');
%[LL1,LH1,HL1,HH1]=dwtper2(X,'bior3.3');
%load lena
X=imread('E:\my\畢業(yè)設(shè)計(jì)\matlab\LENA256.bmp','bmp');
%figure(100);imshow(uint8(X));title('original image');
%利用正交小波函數(shù)對(duì)X進(jìn)行多尺度分解
%[C,S]=wavedec2(double(X),N,Wavelet_Name)是使用指定的小波基‘Wavelet_Name’對(duì)圖像X進(jìn)行N層分解

Wavelet_Name='db1';
Wavelet_Scale=1;

[LL1,LH1,HL1,HH1]=swt2(double(X(:,:,1)),Wavelet_Scale,Wavelet_Name);
%訓(xùn)練LL1頻帶碼本，將碼本大小設(shè)為128，維度為4 %
th=0.05;
nc=128;nd=4;
[cbLL1]=LBG(LL1,nc,nd,th); %LBG.m在L3.3節(jié)中
%訓(xùn)練LH1頻帶碼本，將碼本大小設(shè)為64，維度為4 %
%利用一層小波變換加上向量量化進(jìn)行Lena圖像壓縮%
%輸入一幅Lena圖像%
X=imread('E:\my\畢業(yè)設(shè)計(jì)\matlab\LENA256.bmp','bmp');
%load lena
[M,N]=size(X);
%使用Matlab中Wavelet工具箱的周期性小波變換函數(shù)dwtper2.m%
%并使用雙正交小波bior3.3進(jìn)行一層小波變換%
%[LL1,LH1,HL1,HH1]=dwtper2(X,'bior3.3');
[LL1,LH1,HL1,HH1]=swt2(double(X(:,:,1)),Wavelet_Scale,Wavelet_Name);
wt_X(1:M/2,1:N/2)=LL1;
wt_X(1:M/2,N/2+1:N)=LH1;
wt_X(M/2+1:M,1:N/2)=HL1;
wt_X(M/2+1:M,N/2+1:N)=HH1;
%分別針對(duì)不同的頻帶進(jìn)行向量量化，其中VQ.m置于L3.4節(jié)中%


[R_LL1,bitLL1]=VQ(LL1,cbLL1);

%計(jì)算圖像經(jīng)過(guò)向量量化后的壓縮率CR及失真PSNR%
totalbit=bitLL1     %花費(fèi)的總位數(shù)
%對(duì)高頻的域值處理
%程序L7_2.m: 應(yīng)用小波變換與向量量化做圖像壓縮%
%利用LGB訓(xùn)練一層小波變換后四個(gè)頻帶的碼本%
function Wavelet_Vector_Comprssion()
%X=imread('j:\數(shù)字圖像處理－活用Matlab程序\CH7\lena.bmp','bmp');
%[LL1,LH1,HL1,HH1]=dwtper2(X,'bior3.3');
%load lena
X=imread('E:\my\畢業(yè)設(shè)計(jì)\matlab\LENA256.bmp','bmp');
%figure(100);imshow(uint8(X));title('original image');
%利用正交小波函數(shù)對(duì)X進(jìn)行多尺度分解
%[C,S]=wavedec2(double(X),N,Wavelet_Name)是使用指定的小波基‘Wavelet_Name’對(duì)圖像X進(jìn)行N層分解
Wavelet_Name='db1';
Wavelet_Scale=1;
[C,S]=wavedec2(double(X),Wavelet_Scale,Wavelet_Name);
%appcoef2用來(lái)提取二維信號(hào)小波分解的近似分量
LL1=appcoef2(C,S,Wavelet_Name,Wavelet_Scale);

%detcoef2用來(lái)提取二維信號(hào)小波分解的細(xì)節(jié)分量
%D=detcoef2(O,C,S,N)，O指定細(xì)節(jié)分量的類(lèi)型，
%當(dāng)O=‘h’時(shí)，重構(gòu)水平細(xì)節(jié)分量，當(dāng)O=‘v’時(shí)，重構(gòu)垂直細(xì)節(jié)分量，當(dāng)O=‘d’時(shí)，重構(gòu)細(xì)節(jié)細(xì)節(jié)分量
LH1=detcoef2('v',C,S,1);
HL1=detcoef2('h',C,S,1);
HH1=detcoef2('d',C,S,1);



%訓(xùn)練LL1頻帶碼本，將碼本大小設(shè)為128，維度為4 %
th=0.05;
nc=128;nd=4;
[cbLL1]=LBG(LL1,nc,nd,th); %LBG.m在L3.3節(jié)中
%訓(xùn)練LH1頻帶碼本，將碼本大小設(shè)為64，維度為4 %



%利用一層小波變換加上向量量化進(jìn)行Lena圖像壓縮%
%輸入一幅Lena圖像%
X=imread('E:\my\畢業(yè)設(shè)計(jì)\matlab\LENA256.bmp','bmp');
%load lena
[M,N]=size(X);
%使用Matlab中Wavelet工具箱的周期性小波變換函數(shù)dwtper2.m%
%并使用雙正交小波bior3.3進(jìn)行一層小波變換%
%[LL1,LH1,HL1,HH1]=dwtper2(X,'bior3.3');
[C,S]=wavedec2(double(X),Wavelet_Scale,Wavelet_Name);
LL1=appcoef2(C,S,Wavelet_Name,Wavelet_Scale);
LH1=detcoef2('v',C,S,1);
HL1=detcoef2('h',C,S,1);
HH1=detcoef2('d',C,S,1);

wt_X(1:M/2,1:N/2)=LL1;
wt_X(1:M/2,N/2+1:N)=LH1;
wt_X(M/2+1:M,1:N/2)=HL1;
wt_X(M/2+1:M,N/2+1:N)=HH1;
%分別針對(duì)不同的頻帶進(jìn)行向量量化，其中VQ.m置于L3.4節(jié)中%


[R_LL1,bitLL1]=VQ(LL1,cbLL1);

%計(jì)算圖像經(jīng)過(guò)向量量化后的壓縮率CR及失真PSNR%
totalbit=bitLL1     %花費(fèi)的總位數(shù)
%對(duì)高頻的域值處理
u=mean2(HL1);
sum1=0;
[m,n]=size(HL1)
for i=1:1:m
    for j=1:1:n
   if HL1(i,j)<u
        HL1(i,j)=0;
       sum1=sum1+1; 
   end
    end
end
sum1
u=mean2(LH1);
sum1=0;
[m,n]=size(LH1)
for i=1:1:m
    for j=1:1:n
   if LH1(i,j)<u
       LH1(i,j)=0;
       sum1=sum1+1; 
   end
    end
end
sum1
u=mean2(HH1);
sum1=0;
[m,n]=size(HH1)
for i=1:1:m
    for j=1:1:n
   if HH1(i,j)<u
        HH1(i,j)=0;
       sum1=sum1+1; 
   end
    end
end
sum1
%重構(gòu)還原圖像
Y=iswt2(R_LL1, HL1 ,LH1, HH1,'db1',size(X));
figure(200);imshow(uint8(Y));title('compressed image');
figure(100);imshow(uint8(X));title('original image');

X=double(X);
Y=double(Y);

MSE=(sum(sum((X-Y).^2)))/(M*N)
PSNR=20*log10(255/sqrt(MSE))
CR=M*N*8/totalbit

%figure(200);imshow(uint8(Y));title('compressed image');

%重構(gòu)還原圖像
Y=iswt2(R_LL1, HL1 ,LH1, HH1,'db1');
figure(200);imshow(uint8(Y));title('compressed image');
figure(100);imshow(uint8(X));title('original image');

X=double(X);
Y=double(Y);

MSE=(sum(sum((X-Y).^2)))/(M*N);
PSNR=20*log10(255/sqrt(MSE));
CR=M*N*8/totalbit

%figure(200);imshow(uint8(Y));title('compressed image');