function [final] =lloyd_max(b,initial,ss)

%  lloyd_max() calculates the b-bit quantization of an image
%              according to the Lloyd-Max Quantization Method

%  usage  --->  QUANT_LM = lloyd_max(b,initial,ss)
%  inputs --->	b: 	   number of bit quantization desired
%		initial:   initial image to be quantized
%		ss:	   the bit level of original image 	 
%  output --->	W:	   the final quantized image
%****************************************************************


% First determine the pdf for the image intensity levels (quantization
%   values).  Perform a 2D histogram of the image, and then sum the number
%   of values at each decision level. Normalize the resulting profile.  This
%   will yield the pdf for the entire image.
 
  [M,N]=size(initial);		% M,N: number of rows,columns in initial
  L=2^b;			% L:  number of intensity levels in final
  Ls=2^ss;			% Ls: number of intensity levels in initial
  Rmin=min(min(initial));	% Rmin: minimum intensity level in initial
  Rmax=max(max(initial));	% Rmax: maximum intensity level in initial
 
  A=Rmax-Rmin;  		% A: dynamic range of image in int. levels
  Q=A/(L-1);                      % Q: quantization step size
 

  hst=zeros(1,Ls);
  for m=1:M
   for n=1:N
    hst(initial(m,n)+1)=hst(initial(m,n)+1)+1;
   end;
  end;
 
  px=hst./sum(hst);
 
  n=1:Ls;

 % figure;
 % plot(n,px);
 % title('pdf of Original Image Intensity');

%***********************************************************************
%****** Compute the Lloyd-Max quantization(q) and decision(d) levels ***

 tolerance = 0.01;  %error tolerance in percent for error convergence
 
 
 for k = 1:(L+1);
  d(k) = (k-1)*(A/L);
 end;
 
 ssum = 0.01;
 lastsum = 0.001; 
 
 while ((~(abs(ssum - lastsum)/lastsum < tolerance)) & ~(ssum == 0))
 
% Find quantization levels
   for k = 1:(L);
     limit = 0;
     for k2 = 1:Ls
       if ((((k2-1)) < d(k+1)) & (((k2-1)) >= d(k)))
         limit = [limit,k2];
       end
     end
     if (length(limit) == 1)
       numint = 0;
       denint = 0;
     else
       lolimit = limit(2);
       hilimit = limit(length(limit));
       fx = px.*(n-1);
       numint = ((d(k+1)-d(k))/(2*(hilimit-lolimit)))*((2*sum(fx((lolimit+1):(hilimit-1)))) + fx(lolimit) + fx(hilimit));
       denint = ((d(k+1)-d(k))/(2*(hilimit-lolimit)))*((2*sum(px((lolimit+1):(hilimit-1)))) + px(lolimit) + px(hilimit));
     end
     if (denint == 0) & (numint == 0)
       q(k) = (d(k) + d(k+1))/2;
     elseif (denint == 0)
       error('Integration error - we messed up!');
       %    If the denom integral is zero but the num integral is not, we've
       % integrated incorrectly.  Halt the whole show.
     else
       q(k) = numint / denint;
     end;
   end
         
% Find new decision levels
        
   for k = 2:L
     d(k) = .5*(q(k)+q(k-1));
   end
 
% Check error convergence
 
   ssum = 0;
         
   for k = 1:(L);
     limit = 0;
     for k2 = 1:Ls
       if ((((k2-1)) < d(k+1)) & (((k2-1)) >= d(k)))
         limit = [limit,k2];
        end
     end
     if (length(limit) == 1)
       integral = 0;
     else
       lolimit = limit(2);
       hilimit = limit(length(limit));
       sigmaqsq = (((n-1) - q(k)).^2).*px;
       integral = ((d(k+1)-d(k))/(2*(hilimit-lolimit)))*((2*sum(sigmaqsq((lolimit+1):(hilimit-1)))) + sigmaqsq(lolimit) + sigmaqsq(hilimit));
     end
     ssum = ssum + integral;
   end
   
   if (ssum == 0)
     ;    
   elseif ~(abs(ssum - lastsum)/lastsum < tolerance)
     lastsum = ssum;
   end  
 
 end
     
% Quantize...

 q

 for m=1:M
  for n=1:N
   for k = 1:L
     if ( (initial(m,n) >= d(k)) & (initial(m,n) < d(k+1)))
       final(m,n)=round(q(k));
     elseif ( (d(k) == d(k+1)) & (initial(m,n) == d(k)) )
       final(m,n)=round(d(k));
     end;
   end;
  end;
 end;