quantization.m

使用matlab实现的矢量量化的LBG算法

clc;                         %清除命令
clear;
disp('training the pictures'); 
disp('     start');
disp('   divide picture   ');
t=0;
time=clock;
Q=input('请输入训练集的个数，值的大小在1-13之间:');
for i=1:Q                  %这里可以改（1-13）
    c=int2str(i);            %把整数转换成字符串
    str=strcat(c,'.bmp');    %连接字符串
    picture=imread(str);     %读图象数据
    picture=double(picture);
    [M,N]=size(picture);     %读取矩阵的行数和列数
    m=bitshift(M,-4);         %移位 
    n=bitshift(N,-4);         %移位 
    for j=1:m
        for k=1:n 
            t=t+1;
            b(1:4,1:4,t)=picture((j-1)*4+1:j*4,(k-1)*4+1:k*4);  %扫描图象保存4*4的块
           
        end
    end
end

%现在开始设计码书
disp('  现在开始设计码书   ');
%利用聚类准则
disp('  利用聚类准则   ');


%)求所有矢量的重心，矢量重心用该矢量中所含像素点灰度值的平均值来表
for i=1:t
    meanofvector(i)=mean(mean(b(1:4,1:4,i)));  
    flagofvector(i)=0;                             %初始化，用来表示是否所有的矢量都已经聚类
end
 
count=0;                             %用来统计分类的个数
pcenter=0;                           %记录中心的位置；
countofcode=0; 
while (sum(sum(flagofvector))~=t)    %看是否全为1，以表示是否聚类完成
   count=count+1;
    
%)取一个矢量，求此矢量重心与其它未聚类矢量重心的距离，找出最小距离对应的矢量作为类覆盖的圆心
 
  temp=find(flagofvector==0);    %找出未聚类的矢量，记住其位置
      
   if(isempty(temp))            break;
   else
       for m=2:length(temp)
           temp1(m)=(meanofvector(temp(1))-meanofvector(temp(m)))^2; 
       end 
       [mindistance,postion]=min(temp1);        %找最小数,并且保存其位置
       flagofvector(temp(1))=1; 
       flagofvector(temp(postion))=1;          
       pcenter=temp(postion);
       class(count).number=count;              %标记类
   end
  
 %以这个最小距离的1．02倍(倍数可以根据实际情况改变)作为半径r，
 %形成一个球覆盖，即根据各矢量重心将相应的矢量归于此类，同时记录类中的个数；
 r=1.02*mindistance;
class(count).numsofvector=0;  %统计这个类中矢量的个数
m=0;
flag1=0;                      %临时变量
               
 for i=1:t
   distance=(meanofvector(pcenter)-meanofvector(i))^2;
    
   if(distance<=r) 
       flag=(flagofvector(i)==1);              %看该矢量是否已经被聚类，用其距离，作为重新划分到另外一个类的标准
       if(flag==1)                            
         for j=1:count-1
             if(class(j).distance<distance)    %已经聚类，看距离是否小于当前距离
               x=find(class(j).postion==i) ;                  %找到位置符合的地方
               class(j).postion(x)=0;                          %从该类中删除
               class(j).numsofvector= class(j).numsofvector-1;   %类的矢量书减1
               flag1=1;
              break;
            end
         end         
      end
     
       if(flag1)
          continue;                    %记录如果这个距离小于以前分过类的距离，则依然保持不变
       end
     m=m+1; 
     class(count).postion(m)=i;        %记录属于该类的矢量  
     class(count).distance=distance;           %记录距离，已作为重新划分的标志
     flagofvector(i)=1;
     class(count).numsofvector=class(count).numsofvector+1;
     class(count).vector(1:4,1:4,m)=b(1:4,1:4,i);
  end 
end             %  for  

%求这个类的质心(质心定义与LBG算法中的相同)，以此得到一个码矢量和其对应的矢量；
 countofcode=countofcode+1;
 [m,n,z]=size(class(count).vector);
for j=1:z
   
    if(j==1) vector(1:4,1:4,countofcode)=class(count).vector(1:4,1:4,j);
    else     vector(1:4,1:4,countofcode)=vector(1:4,1:4,countofcode)+class(count).vector(1:4,1:4,j);
    end
end
   vector(1:4,1:4,countofcode)= vector(1:4,1:4,countofcode)/(z);%求质心
end  
% while
time1=clock;
total=time1-time;
sprintf('the total training time is %f secs\n',total(6))

%test image 
disp('now start testing picture');
testpicture=imread('cman.bmp');
subplot(1,2,1);
imshow(testpicture);
title('orogin picture');

[x,y,z]=size(vector);
 for w=1:z                   %保存码书矢量重心
    temp(w)=mean(mean(vector(1:4,1:4,w))); 
 end

[M,N]=size(testpicture);     %读取矩阵的行数和列数
m=bitshift(M,-2);         %移位 
n=bitshift(N,-2);         %移位 
fid=fopen('compress.txt','wt');
  for j=1:m
        for k=1:n 
           p=mean(mean(testpicture((j-1)*4+1:j*4,(k-1)*4+1:k*4))); 
         for w=1:z
           temp2(w)=(p-temp(w))^2; 
         end 
       [mindistance,postion]=min(temp2);        %找最小数,并且保存其位置
       fprintf(fid,'%d',postion);
       fprintf(fid,'%s',' ');
        end
 end
status=fclose(fid);

fid=fopen('compress.txt','rt'); 
A=fscanf(fid,'%d');
disp('now start decode the code:');i=0;
for j=1:m
   for k=1:n 
     i=i+1;
     imagepicture((j-1)*4+1:j*4,(k-1)*4+1:k*4)= vector(1:4,1:4,A(i));
 end
end

status=fclose(fid);
disp('now diaplay testing picture');
subplot(1,2,2);
imshow(imagepicture,[]);
title('the last  picture');
time2=clock;
total1=time2-time1;
sprintf('the total testing time is %f secs\n',abs(total1(6)))
end