
\chapter{Exemple de travaux pratiques}


Dans tout le TD, on considère une image digitale $a(i,j)$ de taille $256 \times 256$ (i.e. $i,j \in \{1,2,...,256\}$),
à 256 niveaux de gris (i.e. $a(i,j) \in \{1,2,...,256\}$). Commencez, à partir d'une des machines du département, par
choisir votre image parmi celles disponibles dans le répertoire $\sim$gousseau/images/ (évitez l'image {\it mouche.tif}
pour des raisons qui apparaîtront plus loin). Pour la lire sous matlab, tapez

\com{a=imread('image.tif');}

\noindent la variable $a$ est alors une matrices de valeur entières. Visualisez a comme une image (les valeurs
correspondent à des niveaux de gris, 1 est le noir et 256 le blanc)

\com{imagesc(a); colormap gray;}

\noindent puis comme une surface

\com{view(-45+180,30);}

\com{surfl(double(a));}

\com{colormap gray; shading interp; axis equal;}


\section{Manipulation de l'histogramme}

On appelle histogramme de $a$ le vecteur $h(x), 1\leq x\leq 256$
$$h(x)=\#\{  (i,j) / a(i,j)=x  \}.$$
On appelle fonction de distribution de l'image le vecteur de m\^eme taille
$$H(x)=\#\{  (i,j) / a(i,j)\leq x  \}.$$
Visualisez les vecteurs h et H correspondant a votre image

\com{[m,n]=size(a);}

\com{a=double(a);}

\com{b=reshape(a,m*n,1);}

\com{h=hist(b,255);}

\com{h=h/(m*n);}

\com{H=cumsum(h);}

\com{plot(h);figure;plot(H);}

\noindent {\it Egalisation d'histogramme:} on modifie les niveaux de gris de l'image de manière à ce que sa nouvelle
fonction de distribution soit le plus proche possible de l'identité. L'algorithme très simple fourni ci-dessous est une
adaptation en discret de la remarque suivante: si une variable aléatoire $X$ à valeurs dans [0,1] a pour fonction de
distribution $F$, alors la variable aléatoire $F(X)$ a pour fonction de distribution l'identité (le vérifier).
Visualisez la nouvelle fonction de distribution, ainsi que le nouvel histogramme. Qu'en pensez-vous ?

\com{H(256)=1;}

\com{b=zeros(m,n);}

\com{b=round((H(a+1)*255));}


%\section{Convergence au sens des distributions}

%On considère le signal discret $s$=(1 1). En utilisant un résul%tat de vos TD précédents, vers
%quoi doit converger la
%suite des convolutions de ce signal avec lui-m\^eme? Effectuez l'expérience
%numériquement à l'aide de la fonction {\it conv} (quelques convolutions
%suffisent).

%\com{s=[1 1];t=s;}

%\com{t=conv(t,s);  \% quelques fois}

%\com{plot(t,'*');}

%\noindent Qu'en pensez-vous ? Convoluez le résultat avec une
%fonction porte (par exemple le signal (1 1 1 1 1)). Qu'observe-t-on ?

\section{Echantillonnage et transformée de Fourier}


{\noindent \bf Attention:} La transformée de Fourier discrète (TFD) se calcule sous Matlab à l'aide des fonctions {\it
fft} (dimension 1) et {\it fft2} (dimension 2). Comme expliqué à la remarque \ref{TFDMatlab} page ~/pageref{TFDMatlab},
cette fonction ne calcule pas les coefficients $\tilde u_k$, mais les coefficients $\hu_k$ qui sont liés aux $\tilde
u_k$ par la formule \ref{uchapeau} page ~\pageref{uchapeau}. Le passage de l'une a l'autre de ces définitions (en
dimension 1 comme en dimension 2) s'effectue au moyen de la fonction {\it fftshift}. D'autre part, Matlab retourne en
fait les coefficients $n*\hu$ en dimension 1, et $m*n*\hu$ en dimension 2 (tapez \com{help fft}). En résumé, les
coefficients $\tilde u_k$ sont obtenus par la formule suivante:

\com{A=(1/m*n)*fftshift(fft2(a))}

\noindent et la transformée inverse est alors

\com{a=m*n*ifft2(fftshift(A))}

\vskip .5cm

\noindent Visualisez la transformée de votre image au moyen de la commande:

\com{imagesc(log(1+abs(A))); colormap(jet); colorbar }

\noindent Répétez l'expérience avec l'image {\it mouche.tif}.


\noindent Sous-échantillonnez l'image $a$ d'un facteur 2 (i.e. gardez un pixel sur 2 dans chaque direction).
Qu'observe-t-on ?

\com{b=zeros(128);}

\com{for k=1:128;for l=1:128; b(k,l)=a(2*k,2*l); end; end;}

\noindent M\^eme expérience avec l'image de la mouche.

\noindent Effectuez la m\^eme expérience en filtrant préalablement l'image. Vous pouvez par exemple multiplier le
spectre de votre fonction par la fonction indicatrice d'un carré centré de taille deux fois plus petite que celle de
l'image. Le résultat est-il meilleur ? pourquoi ? A quoi sont d\^ues les oscillations  visible sur l'image ?

\section{Zoom par 0-padding (prolongement du spectre par des zéros)}

On se propose de zoomer l'image $a$, par exemple d'un facteur deux dans chaque dimension. L'idée est d'ajouter des
zéros à la TFD de l'image de manière à augmenter sa taille tout en préservant son contenu aux faibles fréquences, puis
de reconstruire une image à partir de cette transformée, référez vous au paragraphe \ref{zoom} page ~\pageref{zoom}.
L'algorithme est le suivant :
\begin{itemize}
\item calculez $A$, la TFD de $a$
\item prolongez $A$ en une image $512 \times 512$ par ajout de zéros aux
hautes fréquences (attention: il est préférable d'utiliser la fonction {\it fftshift} avant de prolonger par des
zéros).
\item calculez $b$, TFD inverse de $A$.
\end{itemize}

Qu'observe-t-on ? En particulier aux bords de l'image ? pourquoi ? Effectuez la m\^eme expérience avec une image
comportant une discontinuité nette (par exemple un carré noir sur fond blanc). M\^eme questions.

\section{Translations}

On cherche à effectuer une translation non entière sur une image, c'est à dire à calculer des valeurs en des points
situés ``entre les pixels''. On utilise la proposition \ref{translation}.

Fabriquez le signal $f=\chi_{[60,80]}$, ($\chi$ étant la fonction indicatrice) considéré comme périodique sur
l'intervalle $[0,128]$, translatez le de $1/2$ en multipliant les coefficients de sa TFD par les exponentielle
adéquates. Visualisez le résultat comme une image (en utilisant simplement la commande {\it imagesc}). Qu'observe-t-on?

\noindent Effectuer une nouvelle translation de $1/2$. Retrouve-t-on le signal original?

\section{Phase et module de la TFD}

Nous allons maintenant illustrer l'importance respective de la phase et du module des coefficients de la TFD. Nous
commen\c{c}ons par échanger les phases des TFD de deux images. Choisissez une image $b$ différente de $a$ (de m\^eme
taille). Calculez les TFD $A$ et $B$ respectivement de $a$ et $b$. Remplacez les coefficients de $B$ par leur modules,
puis multipliez les par la phase des coefficients correspondants de A. Utilisez {\it ifft2} et affichez le résultat.
Essayez avec une image de texture pour $a$, puis pour $b$.

\com{A=fft2(a); B=fft2(b);}

\com{B=abs(B).*exp(i*angle(A));}

\com{c=ifft2(B); imagesc(real(c));}

\noindent Pour illustrer le r\^ole du module des coefficients, on effectue l'expérience suivante: partant d'une image
synthétique simple (un carré, un disque...) on change la phase des coefficients de sa TFD de manière aléatoire, puis on
on calcule l'image correspondante. Attention, lors de la modification de la phase il faut prendre garde à respecter les
symétries de la TFD, ainsi qu'a ne pas modifier les coefficient réels $A(0,0), A(\frac N2,\frac N2), A(\frac N2,0),
A(0,\frac N2)$. Une possibilité est de partir d'un bruit blanc et d'utiliser les phases de sa TFD. Qu'observe-t-on?
Effectuez la m\^eme expérience à partir de la texture {\it strawmat.tif}.

\section{Codes Matlab}

\vskip .5 cm \%visualisation de l'histogramme

a=double(a);

[m,n]=size(a);

b=reshape(a,m*n,1);

h=hist(b,256);

h=h/(m*n);

H=cumsum(h);

plot(h);

figure; plot(H);


\vskip .5 cm \%égalisation

b=zeros(m,n);

b=round((H(a+1)*255));

b=double(b);

c=reshape(b,m*n,1);

h=hist(c,255);

h=h/(m*n);

H=cumsum(h);

figure;plot(h);figure;plot(H);



\vskip .5 cm \%visualisation TFD
%remarque: on n'a pas ici divisé par m*n, de sorte que le vecteur b
b=fft2(a);

b=fftshift(b);

imagesc(log(abs(b)+1)); colormap(jet); colorbar;

\vskip .5 cm \%sous-échantillonage brutal de a (d'un facteur 2)

M=floor(m/2); N=floor(n/2);

b=zeros(M,N);

for i=1:M;

for j=1:N;

b(i,j)=a(2*i,2*j);

end

end

imagesc(b);colormap gray;

\vskip .5cm \%idem après filtrage (on suppose ici que m et n sont divisibles par 4)

b=fft2(a);

b=fftshift(b);

c=zeros(m,n);

M4=m/4; N4=n/4; M2=m/2; N2=n/2;

c(M4+1:3*M4,N4+1:3*N4)=b(M4+1:3*M4,N4+1:3*N4);

d=real(ifft2(fftshift(c)));

b=zeros(M2,N2);

for i=1:M2;

for j=1:N2;

b(i,j)=d(2*i,2*j);

end

end

imagesc(b);colormap gray;

vskip .5 cm \%zoom de l'image

\%m et n doivent être pairs

M=m/2; N=n/2;

imagesc(a); colormap gray;

b=fft2(a);

b=fftshift(b);

c=zeros(2*m,2*n);

c(M+1:M+m,N+1:N+n)=4*b(:,:);

\%le facteur 4 est introduit à cause de la définition de la TFD

\%particulière à Matlab (pas de division par la taille de l'image)

c=fftshift(c);

e=ifft2(c);


figure;

imagesc(real(e)); colormap gray;

\vskip .5 cm \%effet de Gibbs par zoom

n=128;

a=zeros(n);

a(n/4:3*n/4,n/4:3*n/4)=255;

imagesc(a);colormap gray;

b=fft2(a);

b=fftshift(b);

c=zeros(2*n);

c(n/2+1:n/2+n,n/2+1:n/2+n)=4*b(:,:);

c=fftshift(c);

e=ifft2(c);

figure;

imagesc(real(e)); colormap gray;

\vskip .5 cm \%translation

taille=128; \% la taille doit être paire

a=zeros(1,taille);

a(taille/4:3*taille/4)=1;

imagesc(a); colormap gray;

\%permet de voir le signal a comme une image ayant toutes ses lignes identiques

b=fftshift(fft(a));

omeg=exp(-2*i*pi/taille);

angle=.5;

T=taille/2;

for n=1:taille;

b(n)=b(n)*omeg $\hat{}$ (((n-1)-T)*angle);

end;

d=ifft(fftshift(b));

figure;

imagesc(real(d));colormap gray;

\vskip .5 cm \%échange de phases

\%b doit être de même taille que a

A=fft2(a);

B=fft2(b);

B=abs(B).*exp(i*angle(A));

c=ifft2(B);

imagesc(real(c));

\vskip .5 cm \%synthèse de texture a partir d'un carre

a=zeros(256);

a(:,:)=256;

a(120:130,120:130)=0;

imagesc(a);colormap gray;

b=randn(256);

B=fft2(b);

A=fft2(a);

C=abs(A).*exp(i*angle(B));

c=ifft2(C);

figure;imagesc(real(c));colormap gray

\% essayez avec d'autres formes...
