Un principe de compression d’image

1. Un principe de compression d’image Cours S.S.I.I., séance 10, décembre 2013, durée : 50 mn, Jean-Paul Stromboni pour les élèves SI3 Après cette séance, vous devriez savoir : • Décrire le principe de la compression JPEG des images numériques • Expliquer ce qu’on entend par composition fréquentielle d’une image • Utiliser la transformée en cosinus discret (DCT) et son inverse • Donner la signification des coefficients DCT • Compresser une image avec des pertes en négligeant les coefficients DCT les plus faibles Ce cours a été rédigé à partir de plusieurs sources, dont le cours de Pierre Nerzic à aller consulter sur : http://perso.univ-rennes1.fr/pierre.nerzic/IN/Cours/ , cf. la page ‘S3P3 – Codages et compression.html’ TD 10 : Application avec SIVP et les scripts Scilab ‘madct.sci’ et ‘idct.sci’ au besoin On appliquera la DCT, puis une quantification et la DCT inverse, à une image de niveaux de gris tirée d’une image couleur à choisir librement, on veillera simplement à découper une image suffisamment petite (par exemple 64x64 pixels) pour alléger les calculs.

2. Schéma de principe de la compression JPEG (tiré de Wikipédia) Division de l’image en blocs de 8x8 pixels appelés ‘macroblocs’ Séparation de chaque bloc en plans Y (luminance), Cr et Cb (chrominances rouge et bleu), ces deux plans étant sous échantillonnés d’un rapport 2 suivant la hauteur et suivant la largeur, Transformation DCT (Discrete Cosine Transform) de chaque bloc : on obtient 8x8 coefficients de Fourier qui définissent la composition fréquentielle du bloc Quantification des coefficients : les plus faibles en valeur absolue sont annulés ou codés sur un nombre de bits plus faible (le pas de quantification est augmenté) : compression avec pertes. Compression des coefficients restants : codage RLE (Run Length Encoding), codage de Huffman ou VLC (Variable Length Coding)

3. //Frequency components of a signal //---------------------------------- // build a sampled at 1000hz containing // pure frequencies at 50 and 70 Hz sample_rate=1000; t = 0:1/sample_rate:0.6; N=size(t,'*'); //number of samples s=sin(2*%pi*50*t)+sin(2*%pi*70*t+%pi/4); d=dct(s); //dct // zero low energy components D(abs(d)<1)= 0; //compression par seuil size(find(d<>0),'*') //only 20 nonzero coefficients out of 600 clf; plot(s,'b'), plot(dct(d,1),'r') //dct(d,1)dct inverse plot2d(f,d) f=[0:N-1]*sample_rate/(2*N); xtitle('coefficients de dct(s)',...'freq (Hz)', 'd') xgrid La DCT (Discrete Cosine Transform) est une variante de la TFD • A partir d’un exemple fourni par Scilab dans l’aide de la fonction ‘dct’ : • noter que s est de taille 601 échantillons, • d= dct(s) est de même taille, c’est un vecteur réel • d(abs(d))<1)=0 annule toutes les composantes du vecteur d inférieures à 1 en valeur absolue • restent 20 composantes non nulles, d’où le taux de compression C=601/20 ~ 30 environ • dct(d,1) calcule la dct inverse du vecteur d modifié, on obtient le vecteur s décompressé

4. Définition de la Discrete Cosine Transform (variante DCT II) • s= [sn, n=0 .. N-1], avec sn=s(n/fe) • D= dct(s)= [Dm, m= 0 .. N-1], avec Dm=D(fm), fm= (m/N)*(fe/2) • pour trouver fm, on note que : • La DCT inverse retrouve s= idct(D)= [sn, n=0 .. N-1], avec sn=s(n/fe) • On voit ci-dessus que sn est décomposé par DCT en une somme de N fonctions périodiques de fréquences différentes, les :

5. Comparaison de DCT et de FFT • s= [sn, n=0 .. N-1], de taille N, avec sn=s(n/fe) • S=fft(s)=[Sk, k=0..N-1], de taille N, avec Sk=s(kfe/N) • D= dct(s)= [Dm, m= 0 .. N-1], avec Dm=D(fm) • S et D ont la même taille que s, soit N • L’exponentielle complexe devient un cosinus, D est un vecteur réel, à la différence des Sk qui sont complexes • Les fréquences des valeurs calculées sont différentes • Pour fft, Sk=S(fk) • Pour dct, Dm=D(fm) • DCT est donc une variante de FFT, qui calcule un spectre réel et non complexe pour des fréquences entre 0 et fe/2 et non pas 0 et fe

6. Si on considère les Dk, k=0.. N-1 de l’exemple précédent, on constate qu’ils sont souvent faibles voire nuls. Pour compresser : On peut supprimer les plus faibles, c’est ce qu’on a fait ici avec : d(abs(d)<1)=0; On peut coder les plus faibles sur moins de bits, c’est-à-dire augmenter le pas de quantification Dans les deux cas, le signal est modifié : compression avec pertes Utilisation de la DCT pour compresser un signal (avec des pertes) size(find(d<>0),'*') //only 20 nonzero coefficients out of 600 // donc C= 600/20=30 f=[0:N-1]*sample_rate/(2*N); xtitle('coefficients de dct(s)','freq (Hz)', 'd') xgrid

7. Utilisation de la DCT pour compresser une image (cet exemple est tiré de l’aide de Scilab : - - > help dct) • Pour l’image A, la DCT 2D est utilisée, et la fréquence se transforme en résolutions, ou fréquences spatiales, verticale et horizontale. • La ligne ‘size(find(d<>0),’*’)’ trouve 165 composantes non nulles sur 1680 • Le taux de compression vaut un peu plus de 10, C= 1680/165

8. Représenter une fréquence spatiale horizontale • i(x)= 0.5+ 0.5*cos(2*%pi*f*x) • x varie de 0 à L • Définition : N pixels entre 0 et L • Période échantillonnage : L/N • Résolution horizontale: fx = N/L • Pixellisation : x= k*L/N, k= 0 … N-1 • ik= 0.5*(1+cos(2*%pi*f*k*L/N) • Normalisation de fe : L/N= 1 • fx=1 pixel par unité de longueur •  x= 0 .. N-1 • ik=0.5*(1+cos(2*%pi*f*x)

9. Utilisation de la DCT dans le principe de compression JPEG • L’image à compresser est découpée en blocs de 8x8 pixels, auxquels on applique la DCT, qui calcule 8x8 coefficients pour chaque bloc selon la formule suivante : • N=8, pixel(x,y), avec x=0..7, et y=0.. 7 est un bloc de 64 pixels • DCT(i,j), i=0..7, j=0..7 est le tableau des 64 coefficients DCT du bloc • C(i) vaut 1 pour i non nul, et sqrt(2) pour i = 0 (de même pour C(j) et j) • Le tableau DCT(i,j), i=0..7, j=0..7 contient le spectre du bloc de pixels, les fréquences spatiales normalisées varient entre 0 et fx/2=0.5 horizontalement et verticalement. • 64 intensités donnent 64 coefficients DCT, taux de compression : C=1 • La DCT inverse reconstitue le bloc de pixels à partir des coefficients DCT(i,j) i=0..7, j=0..7

10. La DCT décompose chaque bloc de pixels en une somme de 64 blocs élémentaires contenant chacun une fréquence spatiale horizontale et une fréquence verticale On voit que la DCT inverse reconstitue le bloc de pixels en faisant une somme pondérée de blocs ou images élémentaires. Une image élémentaire contient une fréquence horizontale et une fréquence verticale. Les coefficients de pondération sont les DCT(i,j), i=0..7, j=0 .. 7. La correspondance entre les indices i et j et la fréquence spatiale normalisée est : fi = i/(2*N), et fj= j/(2N). Origine des fréquences j=7 j=1 j=0 j i=0 • le bloc élémentaire associé au coefficient DCT(i,j), i=0..7, j=0..7 • fréquence horizontale j/16 • fréquence verticale i/16 • Les DCT(i,j) décroissent quand i et j augmentent dans la plupart des cas (les composantes des blocs de pixels ont des fréquences spatiales basses) i=1 DCT(i,j) forts DCT(i,j) faibles i=7 i Tiré du cours de Pierre Nerzic cité page 1

11. Pour compresser un bloc, il faut négliger des coefficients DCT • Les coefficients DCT tendent à s’annuler quand la fréquence augmente, on parcourt le tableau DCT en zigzag pour s’éloigner progressivement de l’origine des fréquences en haut et à gauche : • Codage en zigzag : on parcourt les coefficients DCT(i,j) à (i+j) croissant : • i=j=0, le coefficient DCT(0,0) • i=1, j=0, puis i=0, j=1, les coefficients DCT(1,0) puis DCT(0,1), • i=2, j=0, puis i=1, j=1, puis i=0, j=2, etc … • On note le dernier coefficient non nul, soit K, et on transmet uniquement les K coefficients non nuls : le taux de compression du bloc est alors C=64/K j j=0 j=1 j=7 0,0 i=0 … i=1 i=7 i

12. Pour augmenter le taux de compression, il faut annuler les coefficients les plus faibles, ou augmenter le pas de quantification utilisé pour les coder en binaire j j=0 j=1 • Compression par zonage : on décrit la matrice D des coefficients en ‘zigzag’, et on annule les coefficients à partir d’une distance fixée : • à partir du second, pour i+j >0 (C=64) • à partir du quatrième, pour i+j >1 (C=64/3) • à partir du septième (C=64/6), etc … • on peut procéder en multipliant D par la matrice ci contre pour (i+j) > 2 • Compression par seuil : D(abs(D)<seuil)=0; • Compression par quantification : on crée une matrice Q contenant un pas de quantification Q(i,j) pour chacun des coefficients D(i,j), par exemple en faisant : Q(i,j)= 1+ (i+j)*q • On fait la division entière de D par Q, soit Dq= floor(D./Q); • On comptabilise les K coefficients non nuls restants, et C=64/K • On décompresse en multipliant Dq par Q, puis en appliquant la DCT inverse j=7 0,0 i=0 i=0 i=1 i=1 i=7 i j=0 j=1 i=0 i=1 Q

13. Application à l’aide de Scilab : un exemple sur une image de niveaux de gris // affichage de D dans un plot X=size(D,'c'); Y=size(D,'r'); y=Y:-1:1; x=1:X; xset("colormap", jetcolormap(256)) grayplot(x, y, D, strf="041"); colorbar(min(abs(D)), max(abs(D))) // quantification par zonage [Q,C]=quantzone(sim,6); ssimrec=zeros(X,Y); for i=1:8:Y-1 for j=1:8:X-1 sim=D(i:i+7,j:j+7); sim=sim.*Q; // ici on approxime dctim=dct(sim,1); ssimrec(i:i+7,j:j+7)=dctim; end end // taux de compression C // affichage de l’image imd dans un plot X=size(ssimrec,'c') Y=size(ssimrec,'r') x=X:-1:1; y=1:Y; xset("colormap",graycolormap(256)) Sgrayplot(y,x,imd',strf="041"); colorbar(0,1) // affichage de l’image décompressée scf(2); xset("colormap",graycolormap(256)) Sgrayplot(y,x,ssimrec',strf="041"); colorbar(0,1) // compressionTest.sce exec('quantzone.sci'); // test de compression fichier='cameraman.png'; im=imread(fichier); imshow(im); // on réduit à un plan de niveaux de gris img=rgb2gray(im); // on découpe une image 64x64 pixel ssim=imcrop(img,[96,32,64,64]); imshow(ssim) // format réel (intensité entre 0. et 1.0) imd=im2double(ssim); imshow(imd) // application de ‘dct’ aux blocs 8x8 // extraits de imd X=size(ssim,'c') Y=size(ssim,'r') D=zeros(X,Y); sim=[ ] dctim=[ ] for i=1:8:Y-1 for j=1:8:X-1 sim=imd(i:i+7,j:j+7); dctim=dct(sim); D(i:i+7,j:j+7)=dctim; end end

14. Effets sur l’image ‘cameraman.png’ pour [Q,C]=quantzone(sim, 6), C= 3.047

15. Compression de plage_bitmap.png (180x240), avec un seuil valant 0.1, on supprime plus des trois quarts des coefficients

16. Principe des codages RLE et VLC ( p. ex. code de Huffman) • Ici, la matrice des coefficients de Fourier décrite en zigzag vaut D=[1, 2 , 2, 3, 3, 0, …, 0] avec 58 ‘0’ terminaux • RLE, pour Run Length Encoding, propose de coder ainsi la matrice des 64 coefficients :D= [<1> 1, <2> 2, <3> 3, <58> 0] soient 8 octets au lieu de 64en JPEG, on utilise même un caractère EoB au lieu du <58> 0 terminal • VLC (pour Variable Length Code), compte les occurrences de chaque code nécessaire et module le nombre de bits en fonction de la fréquence d’apparition des codes, par exemple • ici 0 apparaît 58 fois, • 1 une fois, 2 deux fois et 3 trois fois • on pourrait coder 0 sur deux bits, ‘1O’, puis 3 sur trois bits, ‘110’, 2 par ‘1110’ et 1 par ‘1111’ par exemple