%Quantization of a color image into a fixed number of colors
%The colors are chosen from a given image using the vector quantization algorithm (Generalized Lloyd algorithm)
%	Author: Xiaofeng Xu, Polytechnic University  4/21/2002

%the end threshold of iteration
e=0.1;
NN=256;

%Read color image data from a file
img=double(imread('flowers.tif'));
[H,W,color_no]=size(img);
M=H*W;

imshow(uint8(img));

N=1;

%Initiate codebooks
%Caculate the average gray level for each color
%We begin from three colors.
codebooks(1,1,:)=sum(sum(img))/M;

%Cacluate the MSE between current image block and codebook
sumnormsq=0;
for ii=1:H,
	for jj=1:W,
		temp=img(ii,jj,:)-codebooks(1,1,:);
		sumnormsq=sumnormsq+temp(1,1,1)^2+temp(1,1,2)^2+temp(1,1,3)^2;
	end
end
Msave=sumnormsq/M/3;

%1000 is maximum no. of iterations
for LOOP=1:1000,
   %Use the previous codebooks to get current codebooks.   
	for i=1:N,
		c(0+1,i,:)=(1+e)*codebooks(1,i,:);
		c(0+1,N+i,:)=(1-e)*codebooks(1,i,:);
   end
   %
	N=2*N;
	mMSE(1)=Msave;
	i=0;
   for loop=1:100
      %for each training vector compare with each codeword
      for ii=1:H,
			for jj=1:W,
				minnormsq=inf;
				for n=1:N,
	            %Cacluate the MSE between current image block and codebook
					temp=img(ii,jj,:)-c(i+1,n,:);
					normsq=temp(1,1,1)^2+temp(1,1,2)^2+temp(1,1,3)^2;
               if normsq<minnormsq
                  %record the minimum MSE value and corresponding codebook index
                  minnormsq=normsq;
                  %replace a training vector with its corresponding codeword
						Q(ii,jj,:)=c(i+1,n,:);
					end
				end
			end
		end
      
      %caclulate new center vectors
		for n=1:N,
			sumx(1,1,1:3)=0;
			sumn=0;
			for ii=1:H,
				for jj=1:W,
					if Q(ii,jj,:)==c(i+1,n,:),
						sumx=sumx+img(ii,jj,:);
						sumn=sumn+1;
					end
				end
			end
			c(i+1+1,n,:)=sumx/sumn;
		end
      
      i=i+1;
      %Update the MSE between current image block and codebook
      sumnormsq=0;
		for ii=1:H,
			for jj=1:W,
				temp=img(ii,jj,:)-Q(ii,jj,:);
				sumnormsq=sumnormsq+temp(1,1,1)^2+temp(1,1,2)^2+temp(1,1,3)^2;
			end
		end
		mMSE(i+1)=sumnormsq/M/3;
      
      %if satisfy the error condition of iteration, then stop
		if (mMSE(i-1+1)-mMSE(i+1))/mMSE(i-1+1)<=e,break, end
	end
   
   Msave=mMSE(i+1);
	for n=1:N,
   	codebooks(1,n,:)=c(i+1,n,:);
	end
   
   N
   imshow(uint8(Q));
   
   if N==NN,break,end
	if N==2,imwrite(uint8(Q),'c:\LBG002.tif','TIF'),LBGMSE02=3*Msave;end
	if N==4,imwrite(uint8(Q),'c:\LBG004.tif','TIF'),LBGMSE04=3*Msave;end
	if N==8,imwrite(uint8(Q),'c:\LBG008.tif','TIF'),LBGMSE08=3*Msave;end
	if N==16,imwrite(uint8(Q),'c:\LBG016.tif','TIF'),LBGMSE16=3*Msave;end
	if N==32,imwrite(uint8(Q),'c:\LBG032.tif','TIF'),LBGMSE32=3*Msave;end
	if N==64,imwrite(uint8(Q),'c:\LBG064.tif','TIF'),LBGMSE64=3*Msave;end
	if N==128,imwrite(uint8(Q),'c:\LBG128.tif','TIF'),LBGMSE128=3*Msave;end
end

%Draw Distortion Rate Curves
imwrite(uint8(Q),'c:\LBG256.tif','TIF'),LBGMSE256=3*Msave;
rate=log2([2 4 8 16 32 64 128 256]);
figure,plot(rate,[LBGMSE02,LBGMSE04,LBGMSE08,LBGMSE16,LBGMSE32,LBGMSE64,LBGMSE128,LBGMSE256],'bs')
title('rate-distortion curves'),xlabel('bit rate'),ylabel('mean-squared error')
legend('LBG vector quantizer')