encode8_512.m

基于等均值等方差的图象矢量量化Matlab程序

%矢量X的一些特征值计算
figure(1);
sig=imread('lena1.bmp');
sig=double(sig)/255;
[m_sig,n_sig]=size(sig);
siz_word=8;
siz_book=512;
imshow(sig);

load codebook8_512;
%k=64;
%N=512;

%重新排列图像
num=m_sig/siz_word;
ss=siz_word*siz_word;
nn=num*num;
re_sig=[];
for i=1:m_sig
    for j=1:m_sig
        f1=floor(i./siz_word);
        m1=mod(i,siz_word);
        if m1==0
            m1=siz_word;
            f1=f1-1;
        end
        f2=floor(j./siz_word);
        m2=mod(j,siz_word);
        if m2==0
            m2=siz_word;
            f2=f2-1;
        end
        re_sig(num*f1+f2+1,siz_word*(m1-1)+m2)=sig(i,j);
    end 
end
x=re_sig';

%矢量x均值
mx=[];
for j=1:nn
    sum3=0;
    for i=1:k
        sum3=sum3+x(i,j);
    end
    average=sum3/k;
    mx=[mx,average];
end

%矢量x方差
vx=[];
for j=1:nn
    sum4=0;
    for i=1:k
        sum4=sum4+(x(i,j)-mx(j))^2;        
    end
    v2=sqrt(sum4);
    vx=[vx,v2];
end

%搜索初始匹配码字
tic
Sb=[];
D=[];
P=[];
for i=1:nn
    Sa=abs(mx(i)-mj(1));
    p=1;
    for j=2:N
        arg=abs(mx(i)-mj(j));
        if arg<Sa
            Sa=arg;
            p=j;
        end
    end
    Sb=[Sb,Sa];
    dmin=0;
    for m=1:k
        dmin=dmin+(x(m,i)-y(m,p))^2;
    end
    D=[D,dmin];
    P=[P,p];
end
%编码
for i=1:nn
    codenumber(i)=1;
    Amin=mx(1,i)-sqrt(D(i)/k);
    Amax=mx(1,i)+sqrt(D(i)/k);
    Vmin=vx(i)-sqrt(D(i));
    Vmax=vx(i)+sqrt(D(i));    
    T=P;
    for j=1:N
       if P(i)+j<=N            
          if mj(P(i)+j)>=Amax|mj(P(i)+j)<=Amin
              codenumber(i)=T(i);
              break;
          else
              if vj(P(i)+j)>=Vmax|vj(P(i)+j)<=Vmin
                 codenumber(i)=T(i);
                 break;
              else                        
                 sum6=0;
                 for m=1:k
                     sum6=sum6+(x(m,i)-y(m,P(i)+j))^2;
                     if D(i)<sum6
                        codenumber(i)=T(i);
                        break;                         
                     end
                 end
                 if D(i)>sum6
                    D(i)=sum6;                                
                    Amin=mx(1,i)-sqrt(D(i)/k);
                    Amax=mx(1,i)+sqrt(D(i)/k);
                    Vmin=vx(i)-sqrt(D(i));
                    Vmax=vx(i)+sqrt(D(i));                   
                    codenumber(i)=P(i)+j;
                    T(i)=codenumber(i);
                 else
                     break;
                 end                            
              end
          end
       else
           break;
       end
    end
       
    for j=1:N
       if P(i)-j>=1            
          if mj(P(i)-j)>=Amax|mj(P(i)-j)<=Amin
              codenumber(i)=T(i);
              break;
          else
              if vj(P(i)-j)>=Vmax|vj(P(i)-j)<=Vmin
                 codenumber(i)=T(i);
                 break;
              else                        
                 sum6=0;
                 for m=1:k
                     sum6=sum6+(x(m,i)-y(m,P(i)-j))^2;
                     if D(i)<sum6
                        codenumber(i)=T(i);
                        break;                         
                     end
                 end
                 if D(i)>sum6
                    D(i)=sum6;                                
                    Amin=mx(1,i)-sqrt(D(i)/k);
                    Amax=mx(1,i)+sqrt(D(i)/k);
                    Vmin=vx(i)-sqrt(D(i));
                    Vmax=vx(i)+sqrt(D(i));                   
                    codenumber(i)=P(i)-j;
                    T(i)=codenumber(i);
                 else
                     break;
                 end                            
              end
          end
       else
           break;
       end
    end
end
toc
t=toc

%figure(2);
%imshow(sig);