function quant_spectrum = EQA3YM_quantizer( spectrum, QF )
% Quantize DCT spectrum of an image using precomputed quantizer tables

lumi = [  16 11 10 16 24 40 51 61 ...
          12 12 14 19 26 58 60 55 ...
          14 13 16 24 40 57 69 56 ...
          14 17 22 29 51 87 80 62 ...
          18 22 37 56 68 109 103 77 ...
          24 35 55 64 81 104 113 92 ...
          49 64 78 87 103 121 120 101 ...
          72 92 95 98 112 100 103 99] ...
      * QF;
      
chromi = [17 18 24 47 99 99 99 99 ...
          18 21 26 66 99 99 99 99 ...
          24 26 56 99 99 99 99 99 ...
          47 66 99 99 99 99 99 99 ...
          99 99 99 99 99 99 99 99 ...
          99 99 99 99 99 99 99 99 ...
          99 99 99 99 99 99 99 99 ...
          99 99 99 99 99 99 99 99] ...
      * QF;
block_size = 8;

% Vectorize blocks to quantize them together
Y_values = im2col(spectrum(:,:,1), [block_size block_size], 'distinct');
Cb_values = im2col(spectrum(:,:,2), [block_size block_size], 'distinct');
Cr_values = im2col(spectrum(:,:,3), [block_size block_size], 'distinct');

% Doing quantization
Y_quant = bsxfun(@rdivide, Y_values, lumi');
Cb_quant = bsxfun(@rdivide, Cb_values, chromi');
Cr_quant = bsxfun(@rdivide, Cr_values, chromi');

% Reconstruct blocks from vectors
[ X, Y, C ] = size(spectrum);
quant_spectrum = zeros([ X, Y, C ]);
quant_spectrum(:,:,1) = col2im(Y_quant, [block_size block_size], [ X Y ], 'distinct');
quant_spectrum(:,:,2) = col2im(Cb_quant, [block_size block_size], [ X Y ], 'distinct');
quant_spectrum(:,:,3) = col2im(Cr_quant, [block_size block_size], [ X Y ], 'distinct');

% Apply rounding
quant_spectrum = round(quant_spectrum);

end