%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Function: my_dct_compression

function [Compressed_Image] = my_dct_compression(Input_image, Q_matrix, Color_Downsample)

if Color_Downsample == 1
    size_Image=size(Input_image);
    Input_image=double(rgb2ycbcr(Input_image)); %We are coding the image in the (Y,Cb,Cr) color space
    Block_size=size(Q_matrix,1); % Assuming a square Q matrix

    % Allocate memory for results (use different variables to speed up the process a little bit)
    Input_image_coeff=zeros(size_Image);
    Input_image_coeff_quantized=zeros(size_Image);
    Final_Image=zeros(size_Image);
    Compressed_Image=zeros(size_Image);

    % Obtain the DCT transformation matrix
    DCT_matrix = dctmtx(Block_size);

    for i=1:3, %Each channel
        Input_image_coeff(:,:,i)= blkproc(Input_image(:,:,i),[Block_size Block_size],'P1*x*P2',DCT_matrix,DCT_matrix');        
        Input_image_coeff_quantized(:,:,i) = blkproc(Input_image_coeff(:,:,i),[Block_size Block_size],'round(x./P1).*P1',Q_matrix(:,:,i));
        Final_Image(:,:,i) = blkproc(Input_image_coeff_quantized(:,:,i),[Block_size Block_size],'P1*x*P2',DCT_matrix',DCT_matrix);
    end;

else %Downsample color by a factor of 2 in each dimension
    
    size_Y_Image=size(Input_image(:,:,1));
    size_color_Image=size_Y_Image/2;
    Input_image=double(rgb2ycbcr(Input_image)); %We are coding the image in the (Y,Cb,Cr) color space
    Block_size=size(Q_matrix,1); % Assuming a square Q matrix
    
    %Downsample the color channels by a factor 2 in each dimmension
    Input_Color_image=zeros([size_color_Image 2]); %Allocate memory for the Downsampled Color image
    Input_Color_image(:,:,1)=dyaddown(Input_image(:,:,2),1,'m');
    Input_Color_image(:,:,2)=dyaddown(Input_image(:,:,3),1,'m');
    
    Input_Y_image_coeff=zeros(size_Y_Image); %Allocate memory for the Y DCT coefficients
    Input_Y_image_coeff_quantized=zeros(size_Y_Image); %Allocate memory for the Y DCT coefficients    

    Input_Color_image_coeff=zeros([size_color_Image 2]); %Allocate memory for the Color DCT coefficients
    Input_Color_image_coeff_quantized=zeros([size_color_Image 2]); %Allocate memory for the Color DCT coefficients
    Final_Color_Image=zeros([size_color_Image 2]); %Allocate memory for the Final downsampled Color image
    Final_Image=zeros([size_Y_Image 3]);
    Compressed_Image=zeros([size_Y_Image 3]);

    % Obtain the DCT transformation matrix
    DCT_matrix = dctmtx(Block_size);

    %Luminance channel
    i=1;
    Input_Y_image_coeff= blkproc(Input_image(:,:,i),[Block_size Block_size],'P1*x*P2',DCT_matrix,DCT_matrix');
    Input_Y_image_coeff_quantized = blkproc(Input_Y_image_coeff,[Block_size Block_size],'round(x./P1).*P1',Q_matrix(:,:,i));    
    Final_Image(:,:,i) = blkproc(Input_Y_image_coeff_quantized,[Block_size Block_size],'P1*x*P2',DCT_matrix',DCT_matrix);
    
    for i=2:3, %Each Color channel
        Input_Color_image_coeff(:,:,i-1)= blkproc(Input_Color_image(:,:,i-1),[Block_size Block_size],'P1*x*P2',DCT_matrix,DCT_matrix');        
        Input_Color_image_coeff_quantized(:,:,i-1) = blkproc(Input_Color_image_coeff(:,:,i-1),[Block_size Block_size],'round(x./P1).*P1',Q_matrix(:,:,i));
        Final_Color_Image(:,:,i-1) = blkproc(Input_Color_image_coeff_quantized(:,:,i-1),[Block_size Block_size],'P1*x*P2',DCT_matrix',DCT_matrix);
        Final_Image(:,:,i)=imresize(Final_Color_Image(:,:,i-1),2,'bicubic');
    end;
    
end

Compressed_Image=ycbcr2rgb(uint8(Final_Image));