Image quantization in matlab

Question

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).

what is wrong with the quantization method?

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 > 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

edit: This is for example a the result of DWT->quantizing the LH,HL,HH values ->IDWT. the image is filled with squares.