function G = EQA3YM_DCT( g )
% Perform 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;
% Extend normalization matrix for all blocks of the image
norm = repmat(norm, [X/block_size Y/block_size C]);

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