Institut Galile
Tlcom 3
Travaux pratiques de traitement dimages numriques
Deuxime sance
Institut Galile
2010-2011
G. Dauphin & A. Beghdadi
Prambule
Les programmes ncessaires pour faire ce TP sont huff2norm.m, norm2huff.m, frequency.m qui sont des
fonctions ralises par Giuseppe Ridin ([Link]
Les images sont disponibles dans Matlab
A. Codage sans perte : lalgorithme de Huffman
Introduction - Rappels
Le codage de Huffman est un algorithme de compression de donnes sans perte labor par David Albert
Huffman et publi en 1952. Le code est dtermin partir d'une estimation des probabilits d'apparition
des symboles de source, un code court tant associ aux symboles de source les plus frquents. Le
principe du codage de Huffman repose sur la cration d'un arbre binaire, cest--dire un arbre o chaque
nud a deux fils, le fils de gauche est tiquet par un zro et le fils de droite est tiquet par un 1. Chaque
terminaison reprsente un des symboles coder. Le code associ au symbole est une succession de 0 et
de 1, la succession qui correspond au parcours depuis la racine jusqu la terminaison correspondant au
symbole. Une implmentation du code de Huffman est donne par les programmes huff2norm.m et
norm2huff.m
Cet algorithme est en un certain sens optimal cest--dire que la longueur statistiquement
ncessaire pour coder un des symboles est presque gale lentropie
H ( X ) L( X ) H ( X ) 1
1
pi
o H ( X ) p i log 2
(1)
et L( X ) pi C i et pi est la probabilit du symbole i et C i est la longueur
du code associ au symbole i.
1. Application sur des donnes synthtiques
En utilisant la simulation dun ensemble de donnes alatoires ayant une distribution de
probabilits choisie pralablement suivant une fonction dterministe ou de faon alatoire,
montrez que le codage puis le dcodage de ces donnes se fait sans dgradations et que la relation
(1) est respecte.
Instructions Matlab tester pour raliser le programme :
Simulation de donnes alatoires suivant une loi uniforme valeurs dans un alphabet 256 lments :
>>x=ceil(255*rand(1,2000)) ;
On mesure l'histogramme des donnes x en le mettant sous la forme d'un vecteur x(:) puis en appliquant
la fonction Matlab hist.
>>histogrammeObserve=hist(x(:),0:255);
L'axe des abscisses est l'ensemble des lments de l'alphabet et l'axe des ordonnes est le nombre
d'occurrences.
>>figure(1) ; plot(0:255,histogrammeObserve);
Choix alatoire d'une distribution de probabilit :
>>rand(1,256) ; histogrammeDemande=ans/sum(ans) ;
Distribution de probabilit choisie partir d'une fonction :
>>sin(pi*(1:256)/257) ; histogrammeDemande=ans/sum(ans) ;
A partir de ces histogrammes demands on peut construire les fonctions de rpartitions associes que l'on
note FB (b) P ( B b) qui est valeurs dans [0,1]:
>>figure(1) ; plot(0:255,cumsum(histogrammeDemande));
Soit U une variable alatoire qui suit une loi uniforme sur [0,1]. Alors FB1 (U ) est une variable alatoire
qui suit la loi de B. En effet :
P FB1 (U ) b PU FB (b) FB (b)
Il nous suffit donc de trouver une fonction qui inverse la fonction de rpartition, c'est--dire une fonction
qui retrouve l'indice d'un lment d'un vecteur partir d'une valeur approche de sa valeur et ce en
supposant seulement que ce vecteur est bien ordonn. La fonction CptAsyn ralise ceci, t dsigne les
valeurs approches recherches et A dsigne le vecteur bien ordonn. Dans cette longue instruction, ce qui
apparat comme des guillemets est en fait une juxtaposition de deux apostrophes et est interprt lors de
l'excution en Matlab comme une transposition de la matrice, ce qui apparat comme une apostrophe est
en fait une indication de dbut ou de fin de chane de caractres.
>>CptAsyn=inline('reshape(sum(ones(length(t(:).''),1)*(A(:).'')<=t(:)*ones(1,length(A(:))),2),size(t))','t','A');
Si t est un nombre (et non un vecteur) alors le rsultat de cette fonction est le nombre d'lments
infrieurs ou gaux dans A t. A priori A est rang par ordre croissant, aussi le rsultat est l'indice de t
dans A qui donne une valeur infrieure ou gale t. Cette instruction est alors presque synonyme t=2;
A=[0 3]; find(t>=A,1,'last'),
>>CptAsyn(1,[0 3])
>>CptAsyn(4,[0 3])
>>CptAsyn(-1,[0 3])
Si t est un vecteur ou une matrice alors le rsultat est un vecteur ou une matrice de la mme taille que t et
dont les valeurs sont gales l'application de CptAsyn sur chaque valeur de t.
>> CptAsyn([-1; 4; 1],[0 3])
La variable alatoire peut donc tre simul par
>>x=CptAsyn(rand(1,20000),cumsum(histogrammeDemande));
On peut contrler que la variable alatoire simule est bien conforme la distribution de probabilit
demande :
>>histogrammeObserve= hist(x(:),0:255);
>>histogrammeObserve=histogrammeObserve/sum(histogrammeObserve);
>>figure(1) ; plot(0:255,histogrammeObserve,'b',0:255,histogrammeDemande,'r');
La fonction norm2huff fournit le code correspondent aux donnes dans x suppose valeurs entires entre
0 et 255. info est l'ensemble des informations ncessaires pour dcoder la squence. Le code ainsi obtenu
est une succession d'entiers entre 0 et 255, chacun code sur 8 bits.
>>[code,info]=norm2huff(uint8(x));
Longueur du code moyen par symbole
>>length(code)*8/length(x)
Dans la pratique on est oblige aussi de stocker info, mais alors le lien avec l'entropie est un peu modifi.
Pour calculer l'entropie, on peut utiliser la fonction Matlab entropy ou calculer l'histogramme observ hO
et valuer
>> H=sum(hO(hO~=0)/sum(hO).*log2(sum(hO)./hO(hO~=0)));
2. Application des images numriques bruites et filtres
On se propose, dans ce qui suit, d'analyser les effets de quelques traitements sur le signal image en
tudiant l'entropie associe la distribution des niveaux de gris.
2.1. Choisir une image en niveau de gris. Cette image sera note A. Mettre cette image au format double
et valeurs sur l'ensemble 0255 (pour l'affichage avec imshow, il faudra au pralable diviser 255).
a. Dterminez et tracez l'histogramme des niveaux de gris de cette image. En dduire l'entropie
associe cet histogramme.
b. Ajoutez un bruit blanc gaussien l'image A. Soit B l'image ainsi obtenue. Dterminez et tracez
l'histogramme de B. Dduisez l'entropie associe. Comparez l'image A ; que peut-on en dire ?
c. On applique maintenant un filtre passe-bas l'image A. Soit C l'image rsultante. Calculez
l'histogramme et l'entropie associe. Comparez aux deux cas prcdents. Expliquez les rsultats obtenus.
2.2. Dans chacun des trois cas (A,B et C) appliquez le code de Huffman. Comparez les taux de
compression et l'efficacit.
B. Transforme en cosinus discrte sur une image entire
B1. Application de la DCT globale
La transforme en cosinus discrte est utilise dans la compression jpeg parce quelle permet de
concentrer linformation pertinente sur un petit nombre de coefficients.
Choisissez une image entire en niveaux de gris, appliquez la transforme en cosinus discrte. O
se trouvent, au sens de la valeur absolue, les coefficients levs et les coefficients faibles, si on
voulait lire les coefficients du plus levs au plus faible quel parcours parmi les coefficients
faudrait-il faire ? Calculez le PSNR entre limage originale et limage reconstruite en ne
conservant que les coefficients de la DCT dont les coordonnes vrifient i+j<=n. Reprsentez
lvolution du PSNR en fonction du nombre de composantes non-annules.
Instructions Matlab tester pour raliser le programme :
Calcul des coefficients de la dct.
>>imDct=dct2(im) ;
On cherche visualiser lordonnancement des valeurs des coefficients de la DCT. Pour cela on a recours
aux instructions suivantes qui permettent d'afficher une image aprs galisation d'histogramme, mais en
l'occurrence il ne s'agit pas d'une image mais des composantes de la dcomposition en cosinus discrte.
>>[val,ordre]=sort(abs(imDct(:)));
>>imDctOrdre(ordre)=1:prod(size(im));
>>imDctOrdre=reshape(imDctOrdre,size(im));
>>figure(1); imshow(imDctOrdre/max(max(imDctOrdre)));
Annulation des composantes dont la somme des coordonnes est suprieure n
>>[J,I]=meshgrid(1:size(im,2),1:size(im,1));
>>imDctModifie=imDct;
>>imDctModifie(I+J>n)=0;
>>imDeformee=idct2(imDctModifie) ;
B2. Transforme en cosinus discrte locale
Dans la norme de JPEG, la transforme en cosinus discrte est en fait applique par bloc de 8x8. De fait
ce choix diminue de faon importante la complexit des calculs.
Choisissez une image en niveaux de gris dont le nombre de lignes et de le nombre de colonnes
sont chacun des multiples de 8. Comparez visuellement limage obtenue en annulant une certaine
proportion de coefficients de sa DCT lorsque ces coefficients ont t calculs sur lensemble de
limage et lorsquils ont t calculs par bloc de 8x8. Commentez les diffrences.
Instructions Matlab tester pour raliser le programme :
Rajout des lignes et des colonnes une image de faon obtenir une image dont le nombre de lignes et le
nombre de colonnes sont chacun des multiples de 8 :
>>im=[im;zeros(8*ceil(size(im,1)/8)-size(im,1),size(im,2))];
>>im=[im zeros(size(im,1),8*ceil(size(im,2)/8)-size(im,2))];
Application de la dct par bloc de 8x8 :
>>imDct88=blkproc(im,[8 8],@dct2);
Annulation de coefficients par bloc
>>[jBloc,iBloc]=meshgrid(1 :8,1 :8);
>>iBloc=repmat(iBloc,size(im)/8);
>>jBloc=repmat(jBloc,size(im)/8);
>>imDct88Modifiee=imDct88;
>>imDct88Modifiee(iBloc+jBloc<n)=0;
C. Quantifications des coefficients de la DCT
Lapport de la transforme en cosinus discrte nest pas seulement de concentrer ses valeurs sur un plus
petit nombre de coefficients, cest aussi que visuellement limage est peu dforme lorsquon varie un
tout petit peu ces coefficients. Cest cela qui permet de faire de la compression en quantifiant les
coefficients de la DCT.
Pour des images en niveaux de gris, la norme JPEG consiste dabord transformer les niveaux de
gris en valeurs entre -128 et 127, ensuite calculer par bloc de 8x8 les coefficients de la dct, et enfin
quantifier ces coefficients avec un pas de quantification qui dpend de la position du pixel lintrieur du
bloc 8x8. Ces pas de quantifications sont donns par une matrice.
m=[ 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]
Il est possible de modifier le compromis qualit de la compression/mmoire requise en divisant par un
coefficient appel qualite qui varie entre 0 et 1 : 1 pour la meilleur qualit et 0 pour la plus mauvaise. Ces
coefficients sont ensuite quantifis avec lalgorithme de huffman (au pralable on remet les coefficients
entre 0 et 255).
Pour approcher le taux de compression, code est une succession de 0 et de 1 et donc ce vecteur
require en place mmoire un nombre de bits gale sa longueur. info dtient les informations ncessaires
pour retrouver limage de dpart, il sagit essentiellement de la signification des codes mmoires qui sont
cods dans [Link] sous la forme dun vecteur nul sauf en un petit nombre de cases, les valeurs en
ces cases peuvent tre cods avec 8 bits et la position dans ce vecteur avec 16 bits.
Implmentez sur une image en niveau de gris cet algorithme de compression (niveaux de gris
centrs, calcul des coefficients de la dct par blocs de 8x8, quantification des composantes, codage
huffman des composantes). Proposez une autre matrice m et montrez exprimentalement qu
mme taux de compression, lautre matrice donne une image de qualit similaire ou plus
mauvaise. Le coefficient qualite permet justement dajuster lexprimentation de faon obtenir le
mme taux de compression avec deux matrices m diffrentes.
Instructions Matlab tester pour raliser le programme :
Dclaration de la matrice m
>>m(1,1:8)=[ 16 11 10 16 24 40 51 61];
>>m(2,:)=[ 12 12 14 19 26 58 60 55];
>>m(3,:)=[ 14 13 16 24 40 57 69 56];
>>m(4,:)=[ 14 17 22 29 51 87 80 62];
>>m(5,:)=[ 18 22 37 56 68 109 103 77];
>>m(6,:)=[ 24 35 55 64 81 104 113 92];
>>m(7,:)=[ 49 64 78 87 103 121 120 101];
>>m(8,:)=[ 72 92 95 98 112 100 103 99];
>>qualite=1; m=m/qualite;
Quantification
>>imDctQuantifiee=blkproc(imDct,[8 8],'round(x./P1)',m);
Dcodage
>> imDctNouvelle=imDctQuantifiee.*repmat(m,size(imDctQuantifiee)/8);
Taux de compression
>>bpp= (length(find([Link]~=0))*32+length(code))/prod(size(im));
D. Ordonnancement des coefficients de la DCT
Il est possible de gagner en compression en regroupant les coefficients nuls de la DCT. Pour cela on
rordonne les coefficients de la dct suivant un ordre en zigzag, de cette faon les coefficients nuls sont
la fin. Et on remplace les derniers coefficients nuls de chaque bloc par un coefficient de fin de squence
not EOB (End Of Block). En pratique EOB est une valeur gale 1 C max o C max est la composante
la plus leve.
Le vecteur ligne suivant indique lordre de lecture au sein dun bloc, les composantes indiques dsignent
les cases, elles sont numrotes de 1 64 par incrmentation le long des colonnes (9 dsigne la case
positionne sur la premire ligne et la deuxime colonne).
ordre=[1 9 2 3 10 17 25 18 11 4 5 12 19 26 33
41 34 27 20 13 6 7 14 21 28 35 42 49 57 50
43 36 29 22 15 8 16 23 30 37 44 51 58 59 52
45 38 31 24 32 39 46 53 60 61 54 47 40 48 55
62 63 56 64];
Calculez la place mmoire ainsi conomise pour limage choisie. Calculez le rapport de
compression en bpp (bit par pixel).
Instructions Matlab tester pour raliser le programme :
Rangement des coefficients de la DCT en une succession de colonnes de longueur 64.
>>imDctQuantifieeCol=im2col(imDctQuantifiee,[8 8],'distinct');
Nouvel ordonnancement :
>>ordre=[1 9 2 3 10 17 25 18 11 4 5 12 19 26 33];
>>ordre=[ordre 41 34 27 20 13 6 7 14 21 28 35 42 49 57 50];
>>ordre=[ordre 43 36 29 22 15 8 16 23 30 37 44 51 58 59 52];
>>ordre=[ordre 45 38 31 24 32 39 46 53 60 61 54 47 40 48 55];
>>ordre=[ordre 62 63 56 64];
>>imDctQuantifieeColOrdonnee=imDctQuantifieeCol(ordre,:);
Visualisation de lordonnancement :
>>figure(1) ; imshow(abs(imDctQuantifieeCol)/5) ;
>>figure(2) ; imshow(abs(imDctQuantifieeColOrdonnee/5) ;
Supprimer les zros la fin dun vecteur colonne
>>u=[0 1 1 0 1 0 0 0 0 ] .';
>>v=u(1 :find(u~=0,1,'last')) ;
Retrouver une matrice de taille fixe partir dun vecteur et des marqueurs EOB
>>v=[1 7 2 3 3 7 4 8 7] ;
>>EOB=max(v);
>>indiceFin=find(v==EOB)-1 ;
>>indiceDebut=[1 indiceFin(1:end-1)+2];
>>longueur=indiceFin-indiceDebut+1;
>>w=zeros(64,length(indiceFin));
>>for numBloc=1:length(indiceFin)
>> w(1:longueur(numBloc),numBloc)=v(indiceDebut(numBloc):indiceFin(numBloc));
>>end;
Trouver lordre inverse (cet ordreInverse s'applique exactement de la mme faon que ordre):
>>[a,ordreInverse]=sort(ordre);
Recompose rune image partir des blocs 8x8 rangs en une succession de colonnes de 64 lignes :
>>im=col2im(im,[8 8],[256 256],'distinct');
E. Sous-chantillonnage de la chrominance
Dans la norme JPEG, une image couleur est dabord dcompose en chrominance et en luminance, la
chrominance est ensuite sous-chantillonnne (aprs filtrage par un passe-bas) et ensuite chaque
composante colorimtrique est compresse sparment. Il ny a pas besoin de centrer les composantes de
chrominance, elles sont en effet dj en parti centres.
Indiquez le niveau de compression ainsi atteint.
Instructions Matlab tester pour raliser le programme :
Pour convertir en ycbcr, il faut que im soit une image couleur en format doublee et valeurs entre 0 et 1.
>>imycbcr=rgb2ycbcr(im) ;
Pour sous-chantillonner la chrominance
>> imycbcrSsEchY=imycbcr(:,:,1);
>>imycbcrSsEchCb=blkproc(imycbcr(:,:,2),[1 size(imycbcr,2)],@decimate,2);
>>imycbcrSsEchCr=blkproc(imycbcr(:,:,3),[1 size(imycbcr,2)],@decimate,2);
