gmvq.m

GMVQ算法

function GMVQ(image)
%picture for processing with gray picture;strawberry
[height,width]=size(image);
%--------------------------------------------------------------------------
%-----------------------------------------------------STEP_0:initialization
%assume the gaussian mixture component number is m;
m=9;
%-----------------------------COMPONENTARRAY--------------------------------------------
%  mean     |    variance     |    PRIOR_probability    |  data number   |    buffer   |
%---------------------------------------------------------------------------------------
%COM means Componets
COM=zeros(m,5);
for index=1:m
    COM(index,1)=index*180/m;
    COM(index,2)=2;
    COM(index,3)=1/m;
end
TOTALDISTORTION_OLD=0;
TOTALDISTORTION=1000000;
eta=0;
%------------------------------
%   MINROU    |    INDEX      |
%------------------------------
IMAGE_CLUSTER=zeros(height,width,2);
t=1;
%----------------------------------------------------------------------------------------------------------------------------------------------------
while t<50
    COM_OLD=COM;
    TOTALDISTORTION_OLD=TOTALDISTORTION;
    %----------------------------------------------------STEP_1:-----Find Cells
    %To every pixel(for example,image size of(height*width),find the
    %MAXLIKELIHOOD clusters belong to;
    for i=1:height
        for j=1:width
            if image(i,j)~=0
                MINROU=1000;
                INDEX=0;
                for index=1:m
                    rou=COM(index,3)*(1/(sqrt(2*pi)*COM(index,2)+eps))*exp(-0.5*((double(image(i,j))-COM(index,1))/(COM(index,2)+eps)).^2);
                    rou=-log(rou+eps);
                    if rou<MINROU
                        INDEX=index;   
                        MINROU=rou;
                    end
                end
                IMAGE_CLUSTER(i,j,1)=MINROU;
                IMAGE_CLUSTER(i,j,2)=INDEX;
            end
        end
    end
     %--------------------------------------------------STEP_2:total distortion
     TOTALDISTORTION=0;
     for i=1:height
         for j=1:width
             TOTALDISTORTION=TOTALDISTORTION+IMAGE_CLUSTER(i,j,1);
         end
     end
     %-------------------------------------------------STEP_3:stop or continue?
     if abs((TOTALDISTORTION-TOTALDISTORTION_OLD)/TOTALDISTORTION)>eta
         %update PRIOR_PROBABILITY and CENTROID for every class;
         COM(:,:)=0;
         for i=1:height
             for j=1:width  
                 if image(i,j)~=0
                    index=uint8(IMAGE_CLUSTER(i,j,2));
                     COM(index,4)=COM(index,4)+1;
                    COM(index,5)=COM(index,5)+double(image(i,j));
                 end
             end
         end
         for index=1:m
             if COM(index,4)>0
                COM(index,3)=COM(index,4)/(height*width);
                COM(index,1)=COM(index,5)/COM(index,4);
             end
         end
         %using new PRIOR_PROBABILITY and CENTROIDS,calculating the new VARIANCE
         temp=0;
         for index=1:m
             if(COM(index,4)>0)
                 for i=1:height
                     for j=1:width
                         if IMAGE_CLUSTER(i,j,2)==index
                            temp=temp+(double(image(i,j))-COM(index,1)).^2;
                         end
                     end
                 end
                 COM(index,2)=sqrt(temp/COM(index,4));
             end
         end   
     else 
         break;
     end
          t=t+1;
end
t
COM(:,1:3)
%show histogram;
x=1:180;
gaussian=zeros(1,180);
for index=1:m
    if COM(index,1)~=0
        gaussian=gaussian+COM(index,3)*(1/(sqrt(2*pi)*COM(index,2)+eps))*exp(-0.5*((x-COM(index,1))/(COM(index,2)+eps)).^2);
    end
end
figure,plot(x,gaussian,'r');