title: Image quantization in matlab
body: <p>I created a method in Matlab for quantizing an Image.
  But I don't seem to get a good quantization of the image.
  (For example, when I quantize the image after using DWT and then use IDWT, the image is filled with squares).</p>
  
  <p>what is wrong with the quantization method?</p>
  
  <pre><code>function [imQuant, error] = quantizeImage(imOrig, nQuant, nIter)
  [imQuant, error] = quantizeImageHelper(imOrig, nQuant, nIter);
  end
  
  function [imQuant, error] = quantizeImageHelper(imOrig, nQuant, nIter)
  % Helper function for the quantizeImage function.
  
  imOrig = uint8(round(imOrig*255));
  imHist = imhist(imOrig);
  
  % Calculate number of pixels in the image.
  height = size(imOrig, 1);
  width = size(imOrig, 2);
  numOfPixels = height*width;
  % Initialize the error vector.
  error = zeros(1, nIter);
  % Compute the first division of the intensities to segments.
  z = computeInitialZ(imHist, numOfPixels, nQuant);
  % P(z) - the probability of each intensity in the image.
  pz = imHist/numOfPixels;
  % Compute q, the vector of the new intensities of the quantized image.
  q = qFromZ(z, pz, nQuant);
  % Compute the error of the first z and q we calculated.
  error(1) = calculateError(z, q, pz, nQuant);
  % Compute z,q, and error, and keep improving the results until
  % you reach the wanted number of iteration, or you reach convergence
  % of the errors.
  for i=1:nIter-1
      z = zFromQ(q, nQuant);
      q = qFromZ(z, pz, nQuant);
      error(i+1) = calculateError(z, q, pz, nQuant);
      if (error(i) == error(i+1))
          break
      end
  end
  
  % Calculate the new look up table which contains the new intensities
  % of the quantized image.
  lookUpTable = zeros(1,256);
  for i=1:nQuant
      lookUpTable(z(i)+1:z(i+1)+1) = q(i);
  end
  
  % Quantize the image by mapping each pixel to the new corresponding
  % intensity.
  imQuant = zeros(size(imOrig));
  imQuant(:,:) = lookUpTable(imOrig(:, :)+1);
  imQuant = double(imQuant)/255;
  % Crop the error vector to remove reccuring converging values.
  error = error(error ~= 0);
  
  end
  
  function [z] = computeInitialZ(imHist, numOfPixels, nQuant)
  % Compute the initial division of the intensities to segments.
  % Each segment contains approximately the same number of pixels.
  % Input: imHist - the histogram of the original image.
  %        numOfPixels - number of pixels in the image.
  %        nQuant - the number of intensities the output imQuant image will
  %                 have.
  % Output: z - the initial division of the intensities to segments.
  z = zeros(1, nQuant+1);
  z(1) = 0;
  z(nQuant+1) = 255;
  % Use the cumulative histogram for approximate equal division of the
  % segments.
  imHistCumSum = cumsum(imHist);
  pixelsForSegment = round(numOfPixels/nQuant);
  currentPixelCount = pixelsForSegment;
  for i=2:nQuant
      z(i) = find(imHistCumSum &gt; currentPixelCount, 1)-1;
      currentPixelCount = pixelsForSegment*i;
  end
  end
  
  function [z] = zFromQ(q, nQuant)
  % Compute the division of the intensities to segments by using
  % the formula to compute z out of q.
  % Input: q - the new intensities of the quantized image.
  %        nQuant - the number of intensities the output imQuant image will
  %                 have.
  % Output: z - the division of the intensities to segments.
  z = zeros(1, nQuant+1);
  z(nQuant+1) = 255;
  for i=2:nQuant
      z(i) = floor((q(i-1)+q(i))/2);
  end
  end
  
  function [q] = qFromZ(z, pz, nQuant)
  % Compute the new intensities of the quantized image by using the formula
  % the formula to compute q out of z and the intensity probability function.
  % Input: z - the segment division of the intensities.
  %        pz - the intensity probabilty function. 
  %        nQuant - the number of intensities the output imQuant image will
  %                 have.
  % Output: q - the new intensities of the quantized image.
  q = zeros(1, nQuant);
  pzSum = zeros(1, nQuant);
  zpz = zeros(1, nQuant);
  for i=1:nQuant
      pzSum(i) = sum(pz(z(i)+1:z(i+1)+1));
      zpz(i) = ((z(i)+1):(z(i+1)+1))*(pz(z(i)+1:z(i+1)+1));
      q(i) = floor(zpz(i)/pzSum(i));
  end
  end
  
  function [error] = calculateError(z, q, pz, nQuant)
  % Compute the error values from the z and q values.
  % Input: z - the segment division of the intensities.
  %        q - the new intensities of the quantized image.
  %        pz - the intensity probabilty function. 
  %        nQuant - the number of intensities the output imQuant image will
  %                 have.
  % Output: error - the error of the current z and q computation.
  error = 0;
  for i=1:nQuant
     error = error + ((q(i) - (z(i)+1:z(i+1)+1)).^2)*(pz(z(i)+1:z(i+1)+1));
  end
  end
  </code></pre>
  
  <p>edit: This is for example a the result of DWT->quantizing the LH,HL,HH values ->IDWT. the image is filled with squares.</p>
  
  <p><img src="https://i.stack.imgur.com/deQDN.jpg" alt="enter image description here"></p>