% TP Codages JPEG et MPEG-2 - 3SN-M - 2022
                                               
%--------------------------------------------------------------------------
% Fonction de transformee (directe et inverse) en cosinus discrete par blocs
%--------------------------------------------------------------------------
% I_DCT = DCT2DParBlocs(sens,I,methode,taille_bloc)
%
% sortie  : I_DCT = image de la DCT ou IDCT par blocs
% 
% entrees : sens = sens pour la DCT : 'Direct' ou 'Inverse'
%           I = image avant DCT ou IDCT par blocs
%           methode = methode de calcul de la DCT : 'Matlab' ou 'Rapide'
%           taille_bloc = taille des blocs pour la DCT (ici 8x8)
%--------------------------------------------------------------------------

function I_DCT = DCT2DParBlocs(sens,I,methode,taille_bloc)
    
    if methode == "Matlab"
        if sens == "Direct"
            fun = @(block_struct) dct2(block_struct.data);
        elseif sens == "Inverse"
            fun = @(block_struct) idct2(block_struct.data);
        end
    else
        if sens == "Direct"
            fun = @(block_struct) DCT2Rapide(block_struct.data, taille_bloc);
        elseif sens == "Inverse"
            fun = @(block_struct) IDCT2Rapide(block_struct.data, taille_bloc);
        end
    end

    I_DCT = blockproc(I, [taille_bloc taille_bloc], fun);

end

%--------------------------------------------------------------------------
% Fonction de calcul de transformee en cosinus discrete rapide 
% pour un bloc de taille 8x8
%--------------------------------------------------------------------------
% Bloc_DCT2 = DCT2Rapide(Bloc_Origine, taille_bloc)
%
% sortie  : Bloc_DCT2 = DCT du bloc
% 
% entrees : Bloc_Origine = Bloc d'origine
%           taille_bloc = taille des blocs pour la DCT (ici 8x8)
%--------------------------------------------------------------------------
% https://www.nayuki.io/page/fast-discrete-cosine-transform-algorithms
% https://fr.wikipedia.org/wiki/Transform%C3%A9e_en_cosinus_discr%C3%A8te
function Bloc_DCT = DCTRapide(vector)
    
    % on suppose taille_bloc == 8 
          
    C = cos((1:7)' * pi / 16) / 2;
    M1 = [
        C(4), C(4), C(4), C(4)
        C(2), C(6), -C(6), -C(2)
        C(4), -C(4), -C(4), C(4)
        C(6), -C(2), C(2), -C(6)
    ];
    M2 = [
        C(1), C(3), C(5), C(7)
        C(3), -C(7), -C(1), -C(5)
        C(5), -C(1), C(7), C(3)
        C(7), -C(5), C(3), -C(1)
    ];


    V1 = [
        vector(1) + vector(8)
        vector(2) + vector(7)
        vector(3) + vector(6)
        vector(4) + vector(5)
    ];
    V2 = [
        vector(1) - vector(8)
        vector(2) - vector(7)
        vector(3) - vector(6)
        vector(4) - vector(5)
    ];

    X1 = M1 * V1;
    X2 = M2 * V2;

    Bloc_DCT = [
        X1(1)
        X2(1)
        X1(2)
        X2(2)
        X1(3)
        X2(3)
        X1(4)
        X2(4)
    ];

end

function Bloc_DCT2 = DCT2Rapide(I, taille_bloc)
    
    % on suppose bloc == 8 x 8 

    fun1 = @(block_struct) DCTRapide(block_struct.data);
    DCT1 = blockproc(I, [taille_bloc 1], fun1);

    fun2 = @(block_struct) DCTRapide(block_struct.data)';
    DCT2 = blockproc(DCT1, [1 taille_bloc], fun2);

    Bloc_DCT2 = DCT2;

end

%--------------------------------------------------------------------------
% Fonction de calcul de transformee en cosinus discrete inverse rapide
% pour un bloc de taille 8x8
%--------------------------------------------------------------------------
% Bloc_IDCT2 = IDCT2Rapide(Bloc_DCT2,taille_bloc)
%
% sortie  : Bloc_IDCT2 = Bloc reconstruit par DCT inverse
% 
% entrees : Bloc_DCT2 = DCT du bloc 
%           taille_bloc = taille des blocs pour la DCT (ici 8x8)
%--------------------------------------------------------------------------
% https://www.nayuki.io/page/fast-discrete-cosine-transform-algorithms
% https://fr.wikipedia.org/wiki/Transform%C3%A9e_en_cosinus_discr%C3%A8te

function Bloc_IDCT = IDCTRapide(vector)
    
    % on suppose taille_bloc == 8 
          
    C = cos((1:7)' * pi / 16) / 2;
    M1 = [
        C(4), C(4), C(4), C(4)
        C(2), C(6), -C(6), -C(2)
        C(4), -C(4), -C(4), C(4)
        C(6), -C(2), C(2), -C(6)
    ];
    M2 = [
        C(1), C(3), C(5), C(7)
        C(3), -C(7), -C(1), -C(5)
        C(5), -C(1), C(7), C(3)
        C(7), -C(5), C(3), -C(1)
    ];


    V1 = [
        vector(1)
        vector(3)
        vector(5)
        vector(7)
    ];
    V2 = [
        vector(2)
        vector(4)
        vector(6)
        vector(8)
    ];

    X1 = M1' * V1 + M2 * V2;
    X2 = M1' * V1 - M2 * V2;

    Bloc_IDCT = [
        X1(1)
        X1(2)
        X1(3)
        X1(4)
        X2(4)
        X2(3)
        X2(2)
        X2(1)
    ];

end

function Bloc_IDCT2 = IDCT2Rapide(I, taille_bloc)
    
    % on suppose bloc == 8 x 8 

    fun1 = @(block_struct) IDCTRapide(block_struct.data);
    DCT1 = blockproc(I, [taille_bloc 1], fun1);

    fun2 = @(block_struct) IDCTRapide(block_struct.data)';
    DCT2 = blockproc(DCT1, [1 taille_bloc], fun2);

    Bloc_IDCT2 = DCT2;

end