function g = EQA3YM_IDCT( G )
% Perform inverse DCT on the input image

[X, Y, C] = size(G);
block_size = 8;

% Produce normalization matrix of a block
norm = ones(block_size,block_size);
norm(:,1,:) = 1/sqrt(2);
norm(1,:,:) = 1/sqrt(2);
norm(1,1,:) = 1/2;

% figure(2)
% Produce kernel for DCT
kernel = zeros(block_size^2, block_size^2);
[ u, v ] = meshgrid(0:block_size-1, 0:block_size-1);
for x = 1:block_size
    for y = 1:block_size
        % Calculation the the cosine terms of the sum in the equation for a single pixel...
        cos1 = cos( (2*(x-1)+1) * u * pi / (2*block_size) );
        cos2 = cos( (2*(y-1)+1) * v * pi / (2*block_size) );
        % ... and also normalizing it
        prod = norm .* cos1 .* cos2;
        % Vectorize that matrix
        kernel_vector = im2col(prod,[block_size block_size], 'distinct');
        
%       for the visualization of kernels
%        subplot(8,8,(x-1)*block_size+y)
%        imshow( uint8( 1/4 * prod*128+128 ) )
        
        % The variable 'kernel' will hold the cosine terms for each pixel (one row is for one pixel)
        kernel((x-1)*block_size+y,:) = kernel_vector;
    end
end

% Making image "flat" by moving the images of different color channels next to each other
G_flat = reshape(G, X, Y*C);
% Decompose image to blocks and vectorize blocks
G_vectors = im2col(G_flat,[block_size block_size], 'distinct');
% Compute DCT by calculating the sum in the equation and normalizing it furthermore
g_vectors = 2/block_size * kernel * G_vectors;
% Reshape spectrum from vectors to blocks
g_flat = col2im(g_vectors, [block_size block_size], [X Y*C], 'distinct');
% Arrange back images of different channels to a 3D matrix
g = reshape(g_flat, X, Y, C);

end