commit c5fcc7e80fccaf5624a3a229cefa98420fb3687a Author: Laureηt Date: Sun Jun 25 16:32:14 2023 +0200 init diff --git a/Bloc_DCT2.m b/Bloc_DCT2.m new file mode 100644 index 0000000..e69de29 diff --git a/CodageEntropique.m b/CodageEntropique.m new file mode 100644 index 0000000..8b0bfd1 --- /dev/null +++ b/CodageEntropique.m @@ -0,0 +1,23 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Fonction de calcul d'entropie binaire +%-------------------------------------------------------------------------- +% [poids, H] = CodageEntropique(V_coeff) +% +% sorties : poids = poids du vecteur de donnees encode (en ko) +% H = entropie de la matrice (en bits/pixel) +% +% entree : V_coeff = vecteur dont on souhaite calculer l'entropie +%-------------------------------------------------------------------------- + +function [poids, H] = CodageEntropique(V_coeff) + + [count,~] = groupcounts(V_coeff); + freq = count / length(V_coeff); + H = - sum( freq .* log2(freq) ); + + poids = length(V_coeff) * H / 8 / 1024; + +end diff --git a/CompressionJPEG.m b/CompressionJPEG.m new file mode 100644 index 0000000..eaa1160 --- /dev/null +++ b/CompressionJPEG.m @@ -0,0 +1,121 @@ + +% TP CoCodages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Fonction de compression JPEG d'une image +%-------------------------------------------------------------------------- +% [I_Codee,Poids,Compression,nb_coeffs_AC,nb_coeffs_DC] = +% CompressionJPEG(I_Origine,canal,methode,F_Qualite) +% +% sorties : I_Codee = image de DCT quantifiee +% Poids = poids de l'image d'origine en ko pour les differentes +% etapes de la compression +% Compression = taux de compression final +% nb_coeffs_AC = nombre de coefficients AC dans l'image compressee +% nb_coeffs_DC = nombre de coefficients DC dans l'image compressee +% +% entrees : I_Origine = image originale (ou residuelle) +% canal = canal pour le choix de la table de quantification : +% 'Luminance', 'Chrominance' ou 'Residu' +% methode = methode de calcul de la DCT : 'Matlab' ou 'Rapide' +% F_Qualite = facteur de qualite pour la compression +%-------------------------------------------------------------------------- +% Fonctions a coder/utiliser : DCT2DParBlocs.m +% QuantificationDCT.m +% CodageEntropique.m +% CoefficientsACDC.m +%-------------------------------------------------------------------------- + +function [I_Codee, Poids, Compression, nb_coeffs_AC, nb_coeffs_DC] = CompressionJPEG(I_Origine, canal, methode, F_Qualite) + + taille_bloc = 8; + sens = "Direct"; + + % calcul de la DCT + I_DCT = DCT2DParBlocs(sens, I_Origine, methode, taille_bloc); + + % quantification + I_Codee = QuantificationDCT(sens, I_DCT, canal, F_Qualite, taille_bloc); + + % calcul des coeffs AC DC + [Coeff_AC, Coeff_DC] = CoefficientsACDC(I_Codee, taille_bloc); + + % bidule sur les coeffs + nb_coeffs_AC = length(Coeff_AC); + nb_coeffs_DC = length(Coeff_DC); + + [Poids_AC, ~] = CodageEntropique(Coeff_AC); + [Poids_DC, ~] = CodageEntropique(Coeff_DC); + + Poids.H_JPEG = Poids_AC + Poids_DC; + Poids.Origine = numel(I_Origine) / 1024; + Compression = 100 * (1 - Poids.H_JPEG / Poids.Origine); + +end + +%-------------------------------------------------------------------------- +% Fonction de recuperation des coefficients AC/DC de la DCT par blocs +%-------------------------------------------------------------------------- +%[Coeff_AC_Image, Coeff_DC_Image] = CoefficientsACDC(I_Quant, taille_bloc) +% +% sortie : Coeff_AC = vecteur reunissant tous les coefficients AC de +% l'image jusqu'au dernier coefficient non nul +% (taille variable) +% Coeff_DC = vecteur reunissant tous les coefficients DC de +% l'image (taille fixe) +% +% entree : I_Quant = Image de DCT quantifiee +% taille_bloc = taille des blocs pour la DCT (ici 8x8) +%-------------------------------------------------------------------------- +function [Coeff_AC, Coeff_DC] = CoefficientsACDC(I_Quant, taille_bloc) + + fun = @(block_struct) ParcoursBlocZigzag(block_struct.data, taille_bloc*taille_bloc); + zigzag = blockproc(I_Quant, [taille_bloc taille_bloc], fun); + + % on réarange nos blocs + zigzag = reshape(zigzag', taille_bloc*taille_bloc, [])'; + + % on extrait les coeffs + Coeff_DC = zigzag(:, 1); + Coeff_AC = zigzag(:, 2:end); + + % on suppr les zeros + Coeff_AC_nz = []; + for i = 1:size(Coeff_AC, 1) + last_non_zero = find(Coeff_AC(i, :), 1, 'last'); + next_bloc = [Coeff_AC(i, 1:last_non_zero), 1000]; + Coeff_AC_nz = cat(2, Coeff_AC_nz, next_bloc); + end + Coeff_AC = Coeff_AC_nz'; + +end + +%-------------------------------------------------------------------------- +% Fonction de parcours d'un bloc en zigzag pour recuperer les coefficients +% AC/DC de la DCT +%-------------------------------------------------------------------------- +% Vecteur_zigzag = ParcoursBlocZigzag(Bloc_DCT,nb_pixels) +% +% sortie : Vecteur_zigzag = vecteur des coefficients DC/AC ordonnes du bloc +% +% entrée : Bloc_DCT = DCT du bloc courant +% nb_pixels = nombre de pixels dans le bloc +%-------------------------------------------------------------------------- +function Vecteur_zigzag = ParcoursBlocZigzag(Bloc_DCT,nb_pixels) + % Initialisation du vecteur qui contiendra les coefficients + Vecteur_zigzag = zeros(1,nb_pixels); + % Remplissage en partant du debut et de la fin + for k = 1:nb_pixels/2 + n = ceil(sqrt(2*k+1/4)-0.5); + temp = k - n*(n-1)/2; + if (mod(n,2) < 1) + i = temp; + else + i = n + 1 - temp; + end + j = n + 1 - i; + % Positionnement des coefficients dans le vecteur + Vecteur_zigzag(k) = Bloc_DCT(i,j); + Vecteur_zigzag(65-k) = Bloc_DCT(9-j,9-i); + end +end diff --git a/CompressionMPEG.m b/CompressionMPEG.m new file mode 100644 index 0000000..eb6776f --- /dev/null +++ b/CompressionMPEG.m @@ -0,0 +1,50 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Fonction d'encodage MPEG d'une image (+ reference en JPEG) +%-------------------------------------------------------------------------- +% [IRes_Codee,MVdr,Ir_Codee,Poids,Compression] = ... +% CompressionMPEG(Ic_Originale,Ir_Originale,canal,methode,F_Qualite) +% +% sorties : IRes_Codee = image residuelle (DCT quantifiee) +% MVdr = matrice des vecteurs de deplacements relatifs +% Ir_Codee = image de reference (DCT quantifiee) +% Poids = poids de l'image en ko pour les differentes etapes de +% la compression (inclure les mouvements pour MPEG) +% Compression = taux de compression final +% +% entrees : Ic_Originale = image courante originale +% Ir_Originale = image de reference originale +% canal = canal pour le choix de la table de quantification : +% 'Luminance', 'Chrominance' ou 'Residu' +% methode = methode de calcul de la DCT : 'Matlab' ou 'Rapide' +% F_Qualite = facteur de qualite pour la compression +%-------------------------------------------------------------------------- +% Fonctions a coder/utiliser : EstimationMouvement.m +% PredictionImage.m +% CompressionJPEG.m +%-------------------------------------------------------------------------- + +function [IRes_Codee, MVdr, Ir_Codee, Poids, Compression] = ... + CompressionMPEG(Ic_Originale, Ir_Originale, canal, methode, F_Qualite) + + MVdr = EstimationMouvement(Ic_Originale, Ir_Originale); + + Ip = PredictionImage(Ir_Originale, MVdr); + + Ires = Ic_Originale - Ip; + + [Ir_Codee, ~, ~, ~, ~] = CompressionJPEG(Ir_Originale, "Luminance", methode, F_Qualite); + + [IRes_Codee, Poids_IRes_Codee, ~, ~, ~] = CompressionJPEG(Ires, "Residu", methode, F_Qualite); + + [entropie_mvt, ~] = CodageEntropique(MVdr(:)); + + Poids.Origine = numel(Ic_Originale)/1024; + Poids.MVdr = entropie_mvt; + Poids.H_MPEG = Poids_IRes_Codee.H_JPEG + entropie_mvt; + + Compression = 100 * (1 - Poids.H_MPEG / Poids.Origine); + +end diff --git a/DCT2DParBlocs.m b/DCT2DParBlocs.m new file mode 100644 index 0000000..ae2f6cb --- /dev/null +++ b/DCT2DParBlocs.m @@ -0,0 +1,186 @@ + +% 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 \ No newline at end of file diff --git a/DecompressionJPEG.m b/DecompressionJPEG.m new file mode 100644 index 0000000..bfdd32b --- /dev/null +++ b/DecompressionJPEG.m @@ -0,0 +1,29 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Fonction de decompression/reconstruction JPEG d'une image +%-------------------------------------------------------------------------- +% I_Decodee = DecompressionJPEG(I_Quant,canal,F_Qualite,methode) +% +% sorties : I_Decodee = image reconstruite par quantification et DCT inverses +% +% entrees : I_Codee = image de DCT quantifiee +% canal = canal pour le choix de la table de quantification : +% 'Luminance', 'Chrominance' ou 'Residu' +% methode = methode de calcul de la IDCT : 'Matlab' ou 'Rapide' +% F_Qualite = facteur de qualite pour la compression +%-------------------------------------------------------------------------- +% Fonctions a coder/utiliser : QuantificationDCT.m +% DCT2DParBlocs.m +%-------------------------------------------------------------------------- + +function I_Decodee = DecompressionJPEG(I_Codee, canal, methode, F_Qualite) + + taille_bloc = 8; + sens = "Inverse"; + + I_Quant = QuantificationDCT(sens, I_Codee, canal, F_Qualite, taille_bloc); + I_Decodee = DCT2DParBlocs(sens, I_Quant, methode, taille_bloc); + +end \ No newline at end of file diff --git a/DecompressionMPEG.m b/DecompressionMPEG.m new file mode 100644 index 0000000..fc26629 --- /dev/null +++ b/DecompressionMPEG.m @@ -0,0 +1,36 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Fonction d'encodage MPEG d'une image (+ reference en JPEG) +%-------------------------------------------------------------------------- +% [Ic_Decodee,Ir_Decodee] = ... +% DecompressionMPEG(IRes_Codee,Ir_Codee,MVdr,canal,methode,F_Qualite) +% +% sorties : Ic_Decodee = image courante reconstruite +% Ir_Decodee = image de reference reconstruite +% +% entrees : IRes_Codee = image residuelle de DCT quantifiee +% Ir_Codee = image de reference de DCT quantifiee +% MVdr = matrice des vecteurs de deplacements relatifs +% canal = canal pour le choix de la table de quantification : +% 'Luminance', 'Chrominance' ou 'Residu' +% methode = methode de calcul de la DCT : 'Matlab' ou 'Rapide' +% F_Qualite = facteur de qualite pour la compression +%-------------------------------------------------------------------------- +% Fonctions a coder/utiliser : DecompressionJPEG.m +% PredictionImage.m +%-------------------------------------------------------------------------- + +function [Ic_Decodee, Ir_Decodee] = ... + DecompressionMPEG(IRes_Codee, Ir_Codee, MVdr, canal, methode, F_Qualite) + + Ir_Decodee = DecompressionJPEG(Ir_Codee, "Luminance", methode, F_Qualite); + + Ip = PredictionImage(Ir_Decodee, MVdr); + + IdRes = DecompressionJPEG(IRes_Codee, "Residu", methode, F_Qualite); + + Ic_Decodee = Ip + IdRes; + +end diff --git a/Donnees_TP_MPEG-2.mat b/Donnees_TP_MPEG-2.mat new file mode 100644 index 0000000..a6ecea1 Binary files /dev/null and b/Donnees_TP_MPEG-2.mat differ diff --git a/EstimationMouvement.m b/EstimationMouvement.m new file mode 100644 index 0000000..1dd2ef3 --- /dev/null +++ b/EstimationMouvement.m @@ -0,0 +1,113 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Fonction d'estimation du mouvement par "block-matching" +%-------------------------------------------------------------------------- +% MVr = EstimationMouvement(Ic,Ir) +% +% sortie : MVdr = matrice des vecteurs de deplacements relatifs +% +% entrees : Ic = image courante +% Ir = image de reference +%-------------------------------------------------------------------------- + +function MVr = EstimationMouvement(Ic, Ir) + + fun = @(bloc) CSAMacroBloc(bloc.data, bloc.location, Ir); + MVr = blockproc(Ic, [16, 16], fun); + +end + + +%-------------------------------------------------------------------------- +% Fonction de recherche par 'Cross Search Algorithm' : +% - Recherche pour un macro-bloc de l'image courante +%-------------------------------------------------------------------------- +% [MB_IP,Vm] = CSAMacroBloc(MBc, Vp, Iref) +% +% sorties : MB_IP = macro-bloc de prediction dans Iref +% Vm = vecteur de mouvement +% +% entrées : Mbc = macro-bloc dans l'image courante Ic +% Vp = vecteur de prediction (point de depart du MacroBloc) +% Iref = image de reference (qui sera conservee dans le GOP) +%-------------------------------------------------------------------------- + +function Vm = CSAMacroBloc(MB_IC, Vp, IRef) + + %-------------------------------------------------------------------------% + % Initialisation % + %-------------------------------------------------------------------------% + + Vm = zeros(1, 1, 2); + coordx = Vp(1):Vp(1)+15; + coordy = Vp(2):Vp(2)+15; + + %-------------------------------------------------------------------------% + % Vérification du minimum des 9 EQM et déplacement ou arrêt % + % en fonction de la position du minimum parmi les 9 % + %-------------------------------------------------------------------------% + + while true + +% directions = [ +% 0, 0 +% -1, 0 +% 0, -1 +% 1, 0 +% 0, 1 +% ]; % 5 directions + directions = fullfact([3, 3]) - 2; % 9 directions + + fun = @(direction) EQMMacrocBlocVoisin(MB_IC, IRef, size(IRef, 1), size(IRef, 2), coordx, coordy, direction.data); + EQMs = blockproc(directions, [1, 2], fun); + + [~, index] = min(EQMs); + direction_min = directions(index, :); + + if isequal(direction_min, [0, 0]) + break + end + + coordx = coordx + direction_min(1); + coordy = coordy + direction_min(2); + Vm(1, 1, 1) = Vm(1, 1, 1) + direction_min(1); + Vm(1, 1, 2) = Vm(1, 1, 2) + direction_min(2); + + end + +end + +%-------------------------------------------------------------------------- +% Fonction de calcul de l'EQM avec differents voisins +% dans l'image de reference +%-------------------------------------------------------------------------- +% EQM = EQMMacrocBlocVoisin(MB_IC_V,IRef,size_x_Ref,size_y_Ref,coordx,coordy,voisin) +% +% sortie : EQM = erreur quadratique moyenne entre macro-blocs +% +% entrées : MB_IC_V = macro-bloc dans l'image courante (vectorise) +% Ir = Image de reference +% size_x_Ir = nombre de lignes de Ir (pour effets de bords) +% size_y_Ir = nombre de colonnes de Ir (pour effets de bords) +% coordx = les 16 coordonnees du bloc suivant x +% coordy = les 16 coordonnees du bloc suivant y +% voisin = choix du voisin pour decaler le macro-bloc dans Ir +% ('haut', 'gauche', 'centre', 'droite', bas', ...) +%-------------------------------------------------------------------------- + + +function EQM = EQMMacrocBlocVoisin(MB_IC_V, Ir, size_x_Ir, size_y_Ir, coordx, coordy, voisin) + + coordx = coordx + voisin(1); + coordy = coordy + voisin(2); + + if sum((coordy < 1) + (coordx < 1) + (coordy > size_y_Ir) + (coordx > size_x_Ir)) > 0 + EQM = +inf; + else + patch = Ir(coordx, coordy); + EQM = mean((MB_IC_V(:) - patch(:)).^2); + end + +end diff --git a/Exercice_1.m b/Exercice_1.m new file mode 100644 index 0000000..c4e24ff --- /dev/null +++ b/Exercice_1.m @@ -0,0 +1,44 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Exercie_1 : Test du codage entropique +%-------------------------------------------------------------------------- +% Fonction a coder/utiliser : CodageEntropique.m +%-------------------------------------------------------------------------- + +clear +close all +clc +taille_ecran = get(0,'ScreenSize'); +L = taille_ecran(3); +H = taille_ecran(4); +figure('Name','Test du codage entropique',... + 'Position',[0.2*L,0.05*H,0.6*L,0.7*H]); + +%-------------------------------------------------------------------------- + +% Chargement des images de test +I_Ent_1 = load('Donnees_TP_MPEG-2.mat').I_Ent_1; +I_Ent_2 = load('Donnees_TP_MPEG-2.mat').I_Ent_2; +% Calcul de l'entropie avec la fonction Matlab +% (seulement pour une image en niveaux de gris) +H_ref = entropy(I_Ent_1); +% Calcul de l'entropie generalisee pour n'importe quel vecteur +[P,H_v] = CodageEntropique(I_Ent_1(:)); +% Verification de la bonne valeur de l'entropie +if (abs(H_v - H_ref) < 1e-10 && abs(P - 26.9471949267) < 1e-10) + imagesc(I_Ent_1) + colormap gray + axis image off + title({'Le calcul d''entropie est bon :' ... + ['Entropie = ' num2str(H_v,'%.3g') ' bits/pixel (attendue : ' num2str(H_ref,'%.3g') ' bits/pixel)'] ... + ['Poids de l''image = ' num2str(P,'%.3g') ' ko (attendu : ' num2str(26.9) ' ko)']}) +else + imagesc(I_Ent_2) + colormap gray + axis image off + title({'Le codage entropique n''est pas bon :' ... + ['Entropie = ' num2str(H_v,'%.3g') ' bits/pixel (attendue : ' num2str(H_ref,'%.3g') ' bits/pixel)'] ... + ['Poids de l''image = ' num2str(P,'%.3g') ' ko (attendu : ' num2str(26.9) ' ko)']}) +end diff --git a/Exercice_2.m b/Exercice_2.m new file mode 100644 index 0000000..617f5c4 --- /dev/null +++ b/Exercice_2.m @@ -0,0 +1,43 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Exercie_2 : Test de la DCT 2D par blocs +%-------------------------------------------------------------------------- +% Fonction a coder/utiliser : DCT2DParBlocs.m +%-------------------------------------------------------------------------- + +clear +close all +clc +taille_ecran = get(0,'ScreenSize'); +L = taille_ecran(3); +H = taille_ecran(4); +figure('Name','Test de la DCT 2D par blocs',... + 'Position',[0.2*L,0.05*H,0.6*L,0.7*H]); + +%-------------------------------------------------------------------------- + +% Chargement de l'image de test +I_for_DCT = load('Donnees_TP_MPEG-2.mat').I_for_DCT; +% Traitement de la DCT 2D par blocs +taille_bloc = 8; +% Calcul de la DCT 2D par blocs +I_DCT = uint8(DCT2DParBlocs('Direct',I_for_DCT,'Rapide',taille_bloc)); +% Affichage de l'image avant DCT (DCT inverse effectuee) +subplot 121 + imagesc(abs(I_for_DCT)) + colormap gray + axis image off + title('Image avant la DCT (DCT inverse effectuee, valeur absolue)') +% Affichage de l'image apres DCT +subplot 122 + imagesc(I_DCT) + colormap gray + axis image off + if (round(sum(I_DCT(:))/length(I_DCT(:))) == 178) + title('Image apres la DCT : bonne transformation') + else + title('Image apres la DCT : mauvaise transformation') + end + \ No newline at end of file diff --git a/Exercice_3.m b/Exercice_3.m new file mode 100644 index 0000000..164a81e --- /dev/null +++ b/Exercice_3.m @@ -0,0 +1,44 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Exercie_3 : Test de la quantification par blocs +%-------------------------------------------------------------------------- +% Fonction a coder/utiliser : QuantificationDCT.m +%-------------------------------------------------------------------------- + +clear +close all +clc +taille_ecran = get(0,'ScreenSize'); +L = taille_ecran(3); +H = taille_ecran(4); +figure('Name','Test de la quantification par blocs',... + 'Position',[0.2*L,0.05*H,0.6*L,0.7*H]); + +%-------------------------------------------------------------------------- + +% Chargement de l'image de test +I_for_Quant = load('Donnees_TP_MPEG-2.mat').I_for_Quant; +% Facteur de qualite de la quantification +F_Qualite = 1; +% Traitement de la quantification par blocs +taille_bloc = 8; +% Quantification de l'image +I_Quant = uint8(QuantificationDCT('Direct',I_for_Quant,'Luminance',F_Qualite,taille_bloc)); +% Affichage de l'image avant quantification (quantification inverse effectuee) +subplot 121 + imagesc(I_for_Quant) + colormap gray + axis image off + title('Image avant quantification (quantification inverse effectuee)') +% Affichage de l'image apres quantification +subplot 122 + imagesc(I_Quant) + colormap gray + axis image off + if (sum(I_Quant(:)) == 300150) + title('Image apres quantification : bonne quantification') + else + title('Image apres quantification : mauvaise quantification') + end diff --git a/Exercice_4.m b/Exercice_4.m new file mode 100644 index 0000000..fe2680c --- /dev/null +++ b/Exercice_4.m @@ -0,0 +1,49 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Exercice_4 : Compression JPEG et calcul d'entropie +%-------------------------------------------------------------------------- +% Fonction a coder/utiliser : CompressionJPEG.m +%-------------------------------------------------------------------------- + +clear; +close all; +clc; + +taille_ecran = get(0,'ScreenSize'); +L = taille_ecran(3); +H = taille_ecran(4); +figure('Name','Test de la compression JPEG',... + 'Position',[0.2*L,0.05*H,0.6*L,0.7*H]); + +%-------------------------------------------------------------------------- + +% Chargement de l'image de test +I_for_JPEG = load('Donnees_TP_MPEG-2.mat').I_for_JPEG; +% Choix du canal pour la quantification +canal = 'Luminance'; +% Methode de calcul de la DCT 2D par blocs ('Matlab' ou 'Rapide') +methode = 'Matlab'; +% Choix du facteur de qualite +F_Qualite = 30; +% Compression JPEG avec entropie des coefficients AC/DC separes +[I_Codee,Poids,Compression,nb_coeffs_AC,nb_coeffs_DC] = ... + CompressionJPEG(I_for_JPEG,canal,methode,F_Qualite); +% Affichage de l'image avant compression +subplot 121 + imagesc(uint8(I_for_JPEG)) + colormap gray + axis image off + title(['Image d''origine (Poids = ' num2str(Poids.Origine,'%.0f') ' ko)']) +% Affichage de l'image apres DCT et quantification +subplot 122 + imagesc(abs(I_Codee)) + colormap gray + axis image off + title({'Image quantifiee' ... + ['Poids de l''image = ' num2str(Poids.H_JPEG,'%.3g') ' ko (attendu : ' num2str(4.97) ' ko)'] ... + ['Compression = ' num2str(Compression,'%.3g') '% (attendu : ' num2str(92.2) '%)'] ... + ['Nombre de coefficients AC = ' num2str(nb_coeffs_AC) ' (attendu : 12498)'] ... + ['Nombre de coefficients DC = ' num2str(nb_coeffs_DC) ' (attendu : ' num2str(length(I_for_JPEG(:))/64) ')']}) + \ No newline at end of file diff --git a/Exercice_5.m b/Exercice_5.m new file mode 100644 index 0000000..1953f07 --- /dev/null +++ b/Exercice_5.m @@ -0,0 +1,90 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Exercice_5 : Compression/decompression JPEG avec entropie et PSNR +%-------------------------------------------------------------------------- +% Fonctions a coder/utiliser : CompressionJPEG.m +% DecompressionJPEG.m +%-------------------------------------------------------------------------- + +clear; +close all; +clc; + +taille_ecran = get(0,'ScreenSize'); +L = taille_ecran(3); +H = taille_ecran(4); +figure('Name','Test de la compression/decompression JPEG',... + 'Position',[0.1*L,0.05*H,0.8*L,0.7*H]); + +%-------------------------------------------------------------------------- + +% Chargement de l'image de test +I_for_JPEG = load('Donnees_TP_MPEG-2.mat').I_for_JPEG; +% Choix du canal pour la quantification +canal = 'Luminance'; +% Methode de calcul de la DCT 2D par blocs ('Matlab' ou 'Rapide') +methode = 'Rapide'; +% Choix du facteur de qualite (ici vecteur allant de 1% a 97%) +F_Qualite_Max = 97; +F_Qualite = 1:F_Qualite_Max; +% Vecteur pour recuperer le poids de l'image pour chaque facteur de qualite +vecteur_Poids = zeros(1,length(F_Qualite)); +% Vecteur pour recuperer le PSNR de l'image pour chaque facteur de qualite +vecteur_PSNR = zeros(1,length(F_Qualite)); +% Traitement pour chaque facteur de qualite +for f = 1:length(F_Qualite) + % Compression de l'image + [I_Codee,Poids] = CompressionJPEG(I_for_JPEG,canal,methode,f); + % Decompression de l'image + I_Decodee = DecompressionJPEG(I_Codee,canal,methode,f); + % Recuperation du poids des coefficients AC/DC separes + vecteur_Poids(f) = Poids.H_JPEG; + % Calcul et recuperation du PSNR (en uint8 !) + vecteur_PSNR(f) = psnr(uint8(I_Decodee),uint8(I_for_JPEG)); + % Affichage de l'image d'origine (en uint8) + subplot 221 + imagesc(uint8(I_for_JPEG)) + colormap gray + axis image off + title(['Image d''origine (Poids = ' num2str(Poids.Origine,'%.3g') 'ko)']) + % Affichage de l'image reconstruite pour un facteur donne (en uint8) + subplot 122 + imagesc(uint8(I_Decodee)) + colormap gray + axis image off + title({['Image reconstruite pour F_q = ' num2str(f)] ... + ['Poids = ' num2str(Poids.H_JPEG,'%.3g') ' ko'] ... + ['PSNR = ' num2str(vecteur_PSNR(f),'%.3g')]}) + % Courbes du poids et du PSNR en fonction du facteur de qualite + subplot 223 + % Courbe du PSNR + yyaxis left + plot(1:f, vecteur_PSNR(1:f)) + xlim([1 97]) + ylim([10 50]) + xlabel('F_q') + ylabel('PSNR') + % Courbe du poids + yyaxis right + semilogy(1:f, vecteur_Poids(1:f),... + 1:F_Qualite_Max,Poids.Origine*ones(1,F_Qualite_Max),'.') + ylim([0 1.1*Poids.Origine]) + ylabel('Poids (ko)') + legend('PSNR',... + 'Poids de la compression',... + 'Poids d''origine',... + 'Location', 'NorthWest') + grid on + drawnow; + +end +% Affichage de l'image reconstruite pour un facteur donne (en uint8) +subplot 122 + imagesc(uint8(I_Decodee)) + colormap gray + axis image off + title({['Image reconstruite pour F_q = ' num2str(f)] ... + ['Poids = ' num2str(Poids.H_JPEG,'%.3g') ' ko (attendu : 30.7 ko)'] ... + ['PSNR = ' num2str(vecteur_PSNR(f),'%.3g') ' (attendu : 47.4)']}) diff --git a/Exercice_6.m b/Exercice_6.m new file mode 100644 index 0000000..50c4359 --- /dev/null +++ b/Exercice_6.m @@ -0,0 +1,69 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Exercice_6 : Prediction d'image et calcul de residu +%-------------------------------------------------------------------------- +% Fonctions a coder/utiliser : EstimationMouvement.m +% PredictionImage.m +%-------------------------------------------------------------------------- + +clear; +close all; +clc; + +taille_ecran = get(0,'ScreenSize'); +L = taille_ecran(3); +H = taille_ecran(4); +figure('Name','Prediction d''image par Cross Search Algorithm',... + 'Position',[0.1*L,0.05*H,0.8*L,0.7*H]); + +%-------------------------------------------------------------------------- + +% Chargement des images de test (image de reference) +Ir_for_Res = load('Donnees_TP_MPEG-2.mat').Ir_for_Res; +% Affichage de l'image de reference +subplot 231 + imagesc(uint8(Ir_for_Res)) + axis image off + colormap gray + title('Image de reference (Ir)') +% Chargement des images de test (image courante) +Ic_for_Res = load('Donnees_TP_MPEG-2.mat').Ic_for_Res; +% Affichage de l'image courante +subplot 232 + imagesc(uint8(Ic_for_Res)) + axis image off + colormap gray + title('Image courante (Ic)') +% Affichage de l'image de differences +subplot 233 + imagesc(abs(Ir_for_Res-Ic_for_Res)) + axis image off + colormap gray + title('Image de differences (Ir - Ic en valeur absolue)') +% Estimation du mouvement entre les images +MVr = EstimationMouvement(Ic_for_Res, Ir_for_Res); +% Affichage du flux optique +subplot 234 + imagesc(ColorationFluxOptique(MVr(:,:,1),MVr(:,:,2))) + axis image off + title({'Mouvements relatifs des blocs' ... + ['Somme des mouvements : ' num2str(sum(abs(MVr(:)))) ' (attendue : 5553)']}) +% Prediction de l'image courante avec l'image de reference +Ip = PredictionImage(Ir_for_Res, MVr); +% Affichage de l'image predictive +subplot 235 + imagesc(uint8(Ip)) + axis image off + colormap gray + title({'Image predictive (Ip)' ... + ['Niveau de gris moyen : ' num2str(mean(uint8(Ip(:))),'%.4g') ' (attendu : 82.48)']}) +% Calcul de l'image residuelle +Ires = Ic_for_Res - Ip; +% Affichage de l'image residuelle +subplot 236 + imagesc(abs(Ires)) + colormap gray + axis image off + title('Image residuelle (Ip - Ic en valeur absolue)') diff --git a/Exercice_7.m b/Exercice_7.m new file mode 100644 index 0000000..570e2dc --- /dev/null +++ b/Exercice_7.m @@ -0,0 +1,57 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Exercice_7 : Compression MPEG et calcul d'entropie +%-------------------------------------------------------------------------- +% Fonctions a coder / utiliser : CompressionJPEG.m +% CompressionMPEG.m +%-------------------------------------------------------------------------- + +clear; +close all; +clc; + +taille_ecran = get(0,'ScreenSize'); +L = taille_ecran(3); +H = taille_ecran(4); +figure('Name','Test de la compression MPEG',... + 'Position',[0.1*L,0.05*H,0.8*L,0.7*H]); + +%-------------------------------------------------------------------------- + +% Chargement des images de test +Ir_for_MPEG = load('Donnees_TP_MPEG-2.mat').Ir_for_MPEG; +Ic_for_MPEG = load('Donnees_TP_MPEG-2.mat').Ic_for_MPEG; +% Choix du canal pour la quantification +canal = 'Luminance'; +% Methode de calcul de la DCT 2D par blocs ('Matlab' ou 'Rapide') +methode = 'Matlab'; +% Choix du facteur de qualite (ici vecteur allant de 1% a 97%) +F_Qualite = 30; +% Compression JPEG avec entropie des coefficients AC/DC separes +[Ic_Codee,PCou] = CompressionJPEG(Ic_for_MPEG,canal,methode,F_Qualite); +[IRes_Codee,MVdr,Ir_Codee,PRes,Compression] = ... + CompressionMPEG(Ic_for_MPEG,Ir_for_MPEG,canal,methode,F_Qualite); +% Affichage de l'image avant compression +subplot 211 + imagesc(Ic_for_MPEG) + colormap gray + axis image off + title({'Image courante d''origine' ... + ['Poids = ' num2str(PRes.Origine,'%.3g') ' ko (attendu : 413 ko)']}) +% Affichage de l'image courante apres DCT et quantification +subplot 223 + imagesc(abs(Ic_Codee)) + colormap gray + axis image off + title({'Image courante quantifiee' ... + ['Poids = ' num2str(PCou.H_JPEG,'%.3g') ' ko (attendu : 13.3 ko)']}) +% Affichage de l'image residuelle apres DCT et quantification +subplot 224 + imagesc(abs(IRes_Codee)) + colormap gray + axis image off + title({'Image residuelle quantifiee de l''image courante' ... + ['Poids total = ' num2str(PRes.H_MPEG,'%.3g') ' ko (attendu : 5.96 ko)'] ... + ['Poids des mouvements = ' num2str(PRes.MVdr,'%.3g') ' ko (attendu : 0.996 ko)']}) diff --git a/Exercice_8.m b/Exercice_8.m new file mode 100644 index 0000000..f745a80 --- /dev/null +++ b/Exercice_8.m @@ -0,0 +1,99 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Exercice_8 : Codage/decodage MPEG avec entropie et PSNR +%-------------------------------------------------------------------------- +% Fonctions a coder / utiliser : CompressionJPEG.m +% CompressionMPEG.m +% DecompressionJPEG.m +% DecompressionMPEG.m +%-------------------------------------------------------------------------- + +clear; +close all; +clc; + +delete("saves/exo8.gif") + +taille_ecran = get(0,'ScreenSize'); +L = taille_ecran(3); +H = taille_ecran(4); +figure('Name','Test de la compression/decompression MPEG',... + 'Position',[0.1*L,0.05*H,0.8*L,0.7*H]); + +%-------------------------------------------------------------------------- + +% Creation de l'image residuelle +Ir_for_MPEG = load('Donnees_TP_MPEG-2.mat').Ir_for_MPEG; +Ic_for_MPEG = load('Donnees_TP_MPEG-2.mat').Ic_for_MPEG; +% Choix du canal pour la quantification +canal = 'Luminance'; +% Methode de calcul de la DCT 2D par blocs ('Matlab' ou 'Rapide') +methode = 'Matlab'; +% Choix du facteur de qualite (ici vecteur allant de 1% a 97%) +F_Qualite_Max = 97; +F_Qualite = 1:F_Qualite_Max; +% Vecteur pour recuperer le poids de l'image pour chaque facteur de qualite +vecteur_Poids = zeros(2,length(F_Qualite)); +% Vecteur pour recuperer le PSNR de l'image pour chaque facteur de qualite +vecteur_PSNR = zeros(2,length(F_Qualite)); +% Traitement pour chaque facteur de qualite +for f = 1:length(F_Qualite) + % Compression + [Ic_Codee,PCou] = CompressionJPEG(Ic_for_MPEG,canal,methode,f); + [IRes_Codee,MVdr,Ir_Codee,PRes,Compression] = ... + CompressionMPEG(Ic_for_MPEG,Ir_for_MPEG,canal,methode,f); + % Reconstruction + Ic_Decodee_JPEG = DecompressionJPEG(Ic_Codee,canal,methode,f); + [Ic_Decodee_MPEG,Ir_Decodee] = ... + DecompressionMPEG(IRes_Codee,Ir_Codee,MVdr,canal,methode,f); + % Recuperation du poids des coefficients AC/DC separes + vecteur_Poids(1,f) = PCou.H_JPEG; + vecteur_Poids(2,f) = PRes.H_MPEG; + % Calcul et recuperation du PSNR + vecteur_PSNR(1,f) = psnr(uint8(Ic_Decodee_JPEG),uint8(Ic_for_MPEG)); + vecteur_PSNR(2,f) = psnr(uint8(Ic_Decodee_MPEG),uint8(Ic_for_MPEG)); + % Affichage de l'image d'origine + subplot 222 + imagesc(uint8(Ic_Decodee_JPEG)) + colormap gray + axis image off + title({['Image reconstruite (JPEG) pour F_q = ' num2str(f)] ... + ['(Poids = ' num2str(PCou.H_JPEG,'%.3g') ' ko ' ... + 'et PSNR = ' num2str(vecteur_PSNR(1,f),'%.3g') ')']}) + % Affichage de l'image reconstruite pour un facteur donne + subplot 224 + imagesc(uint8(Ic_Decodee_MPEG)) + colormap gray + axis image off + title({['Image reconstruite (MPEG) pour F_q = ' num2str(f)] ... + ['(Poids = ' num2str(PRes.H_MPEG,'%.3g') ' ko ' ... + 'et PSNR = ' num2str(vecteur_PSNR(2,f),'%.3g') ')']}) + % Courbes du poids et du PSNR en fonction du facteur de qualite + subplot 121 + % Courbe du PSNR + yyaxis left + plot(1:f, vecteur_PSNR(:,1:f)) + xlim([1 97]) + ylim([15 55]) + xlabel('F_q') + ylabel('PSNR') + % Courbe du poids + yyaxis right + semilogy(1:f, vecteur_Poids(:,1:f),... + 1:F_Qualite_Max,PCou.Origine*ones(1,F_Qualite_Max),'.') + ylim([0 1.1*PCou.Origine]) + ylabel('Poids (ko)') + legend('PSNR avec compression JPEG',... + 'PSNR avec compression MPEG',... + 'Poids de la compression JPEG',... + 'Poids de la compression MPEG',... + 'Poids d''origine',... + 'Location', 'NorthWest') + grid on + drawnow; + + export_fig(gcf, "saves/exo8.gif", '-png', '-painters', '-m2', '-append'); + +end \ No newline at end of file diff --git a/PredictionImage.m b/PredictionImage.m new file mode 100644 index 0000000..ac5c6cf --- /dev/null +++ b/PredictionImage.m @@ -0,0 +1,27 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Prediction de l'image courante avec l'image de reference et le mouvement +%-------------------------------------------------------------------------- +% Ip = PredictionImage(Ir,MVr) +% +% sortie : Ip = image predictive +% +% entrees : Ir = image de reference +% MVr = matrice des vecteurs de déplacements relatifs +%-------------------------------------------------------------------------- + +function Ip = PredictionImage(Ir,MVr) + + + Ip = zeros(size(Ir)); + for i = 1:size(MVr, 1) + for j = 1:size(MVr, 2) + + Ip( (i-1)*16+1 : i*16 , (j-1)*16+1 : j*16 ) = Ir( (i-1)*16+1+MVr(i,j,1) : i*16+MVr(i,j,1) , (j-1)*16+1+MVr(i,j,2) : j*16+MVr(i,j,2) ); + + end + end + +end diff --git a/QuantificationDCT.m b/QuantificationDCT.m new file mode 100644 index 0000000..034bd39 --- /dev/null +++ b/QuantificationDCT.m @@ -0,0 +1,83 @@ + +% TP Codages JPEG et MPEG-2 - 3SN-M - 2022 + +%-------------------------------------------------------------------------- +% Fonction de quantification (directe et inverse) des coefficients DCT +%-------------------------------------------------------------------------- +% I_Quant = QuantificationDCT(sens,I_DCT,canal,F_Qualite,taille_bloc) +% +% sorties : I_Quant = image quantifiee +% +% entrees : sens = sens de la quantification : 'Direct' ou 'Inverse' +% I_DCT = image avant la quantification +% canal = canal pour le choix de la table de quantification : +% 'Luminance', 'Chrominance' ou 'Residu' +% F_Qualite = facteur de qualite pour la compression +% taille_bloc = taille des blocs pour la DCT (ici 8x8) +%-------------------------------------------------------------------------- + +function I_Quant = QuantificationDCT(sens, I, canal, F_Qualite, taille_bloc) + + T = ChoixTableQuantification(canal, F_Qualite); + + if sens == "Direct" + fun = @(block_struct) round( block_struct.data ./ T ); + elseif sens == "Inverse" + fun = @(block_struct) round( block_struct.data .* T ); + end + + I_Quant = blockproc(I, [taille_bloc taille_bloc], fun); + +end + +%-------------------------------------------------------------------------- +% Fonction de selection d'une table de quatification JPEG en fonction du +% type de canal et du facteur de qualite +%-------------------------------------------------------------------------- +% T_Quant = ChoixTableQuantification(Canal, F_Qualite) +% +% sortie : T_Quant = Table de quantification JPEG +% +% entrees : canal = 'Luminance' ou 'Chrominance' (tables differentes) +% F_Qualite = facteur de qualite (entre 0 et 97) +%-------------------------------------------------------------------------- + +function T_Quant = ChoixTableQuantification(canal,F_Qualite) + + if canal == "Luminance" + T50 = [ + 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 + ]; + elseif canal == "Chrominance" + T50 = [ + 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 + ]; + elseif canal == "Residu" + T50 = ones(8, 8) * 16; + end + + if F_Qualite > 50 + T_Quant = (100 - F_Qualite) * T50 / 50; + elseif F_Qualite < 50 + T_Quant = 50 * T50 / F_Qualite; + else + T_Quant = T50; + end + + T_Quant = round(T_Quant); + +end \ No newline at end of file