clear all; 
data=imread('Cameraman.bmp');      %调入原始图像 
data=double(data)/255;                     %归一化
[m,n]=size(data);                          %求出图像的行数和列数 
figure(1)
subplot(1,2,1);
imshow(data);                              %显示原始图像 
title('原始图像')
subplot(1,2,2);
imhist(data);
title('直方图')
siz_word=4;                                 %设置码字的大小 
siz_book=512;                               %设置码书的大小 
data1=zeros(m*n,1);                          %拉成列向量
for i=1:m
    for j=1:n
        data1((i-1)*n+j)=data(i,j);
    end
end
M1=floor(m*n/siz_word);                      %计算分成M1段
r=mod(m*n,siz_word);
if r>0
    M1=M1+1;
end
data2=zeros(M1,siz_word);                     %分成M1段
l=1;
A=zeros(siz_word,1);
r=1;
for i=1:m*n
    A(r)=data1(i);
    if r==siz_word
        data2(l,:)=A;
        l=l+1;
        r=1;
    else
       r=r+1;
    end
end
code_book=zeros(siz_book,siz_word);


%LBG算法开始
%初始化码书
l=1;
r=1;
A=zeros(siz_word,1);
for i=1:siz_book*siz_word        %前siz_book个向量作为初始码书
    A(r)=data1(i);
    if r==siz_word
        code_book(l,:)=A;
        l=l+1;
        r=1;
    else
        r=r+1;
    end
end
MIU=zeros(M1,siz_book);
 
t=1;
while t==1
   for i=1:M1                        %分类
       B=zeros(siz_word,1);
       B=data2(i,:);
       A=zeros(siz_word,1);
       A=code_book(1,:);
       tep=0.0;
       for l=1:siz_word
           tep=tep+(A(l)-B(l))^2;
       end
       r=1;
       for j=2:siz_book
           A=code_book(j,:);
           temp=sum((A-B).^2);
           if temp<tep
               r=j;
               tep=temp;
           end
       end
       MIU(i,r)=1.0;           %获胜输出的码字（第i个向量 属于r个码字范围）
   end
   t=0;                        %根据MIU，计算码字
   code_book1=zeros(siz_book,siz_word);
   for j=1:siz_book
       for l=1:siz_word
           tep=0.0;
           for i=1:M1
           code_book1(j,l)=code_book1(j,l)+MIU(i,j)*data2(i,l);
           tep=tep+MIU(i,j);
           end
           if tep>0
              code_book1(j,l)=code_book1(j,l)/tep;
           else
              code_book1(j,l)=0.0;
           end
       end
   end
   tep=0.0;
   for j=1:siz_book
       for l=1:siz_word
           tep=tep+(code_book1(j,l)-code_book(j,l))^2;
       end
   end
   if tep/siz_book<0.000001       %判断迭代次数是否满足
       t=0;
   end
   code_book=code_book1;
end
%编码后图像恢复过程
data3=zeros(M1,siz_word);         %寻找每个输入矢量的获胜神经元作为代表
for i=1:M1
    for j=1:siz_book
      if MIU(i,j)==1
          t=j;
      end
    end
    data3(i,:)=code_book(t,:);
end                 
data5=zeros(m,n);
for i=1:m                        %还原成初始输入矩阵形状
    for j=1:n
        tep=(i-1)*n+j;
        i1=floor(tep/siz_word);
        if i1==0
            i1=1;
        end
        j1=mod(tep,siz_word);
        if j1==0
            j1=siz_word;
        end
        data5(i,j)=floor(data3(i1,j1)*255);
    end
end                          
figure(2)
imshow(uint8(data5));  %显示恢复图像
title('矢量量化编码后恢复的图像')
psnr1 = psnr(data5,data)
ssim1 = ssim(data5,data)