%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Feuille de TP n° 3 pour le cours ``Analyse hilbertienne et analyse de Fourier'',%
% de Jean-Michel Morel et Sylvie Fabre.                                           %
% Année scolaire 2009-2010.                                                       %
%                                                                                 %
% 2009/10/16                                                                      %
% Auteur : Bruno Galerne (galerne@cmla.ens-cachan.fr)                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\documentclass[12pt,a4paper]{article}
\usepackage[french]{babel}
\usepackage[T1]{fontenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{lastpage}
\usepackage{times}
\usepackage{stmaryrd}
\usepackage{amsthm}

% Mise en page
\setlength{\textheight}{24.2cm}
\setlength{\textwidth}{16.0cm}
\setlength{\oddsidemargin}{0.0cm}
\setlength{\evensidemargin}{0.0cm}
% \setlength{\parindent}{0.0cm} %
\setlength{\topmargin}{-1cm} %%

% Numerotation des pages / la derniere
\makeatletter
\renewcommand{\@evenfoot}%
        {\hfill{}\thepage \slash \pageref{LastPage}}
\renewcommand{\@oddfoot}{\@evenfoot}
\makeatother

% % Suppression de l'indentation
% \setlength{\parindent}{0pt}

\newcommand{\R}{\mathbb{R}}
\newcommand{\N}{\mathbb{N}}
\renewcommand{\epsilon}{\varepsilon}
\renewcommand{\geq}{\geqslant}
\renewcommand{\leq}{\leqslant}

\def\ds{\displaystyle }


% Définition des theoremes
\newtheorem{theoreme}{Théorèm}
\newtheorem{definition}[theoreme]{Définition}
\newtheorem{proposition}[theoreme]{Proposition}



% Définition de l'environnement question :
\newcounter{questcount} % Compteur des quetions.
\newlength{\longsymb} % Longueur du symbole à gauche des questions.
\settowidth\longsymb{$\blacktriangleright$~}
\newlength{\parshift} % Longueur du shift dans l'indentation des paragraphes des questions.
\setlength{\parshift}{0mm}
\addtolength{\parshift}{\parindent}
\addtolength{\parshift}{\longsymb}
\newenvironment{question}{%
\addtolength{\parindent}{-\parshift} % décalage à gauche
\stepcounter{questcount} % incrémentation du compteur de question
\par % fin du paragraphe précédent
$\blacktriangleright$~\textbf{Question~\thequestcount~:}
}{%
\addtolength{\parindent}{\parshift}
\par % fin du paragraphe courant
}



\begin{document}

{\raggedright Analyse hilbertienne et analyse de Fourier}

{\raggedright ENS de Cachan}

{\raggedright Département de Mathématiques\hfill{}Vendredi 4 décembre 2009}

\vspace{0.5cm}

{\centering \LARGE{\textbf{Feuille de TP numéro 3}}\\}

\vspace{0.5cm}


\section{Consignes et mise en route}


\subsection{Rapport de TP}

\begin{description}
\item[- Consignes pour le rapport de TP :] cette feuille de TP comporte des questions auxquelles il est demandé de répondre dans un rapport de TP.
\item[- Barème :] la moyenne des notes obtenues pour les rapports des quatre TP contribuera à la note de partiel à hauteur de 6 points sur 20.
\item[- Format du rapport :] les rapports de TP sont à rendre sous forme de fichier .pdf rédigés avec \LaTeX.
\item[- Remise du rapport :] le rapport de ce TP doit être \underline{impérativement} rendu avant le \textbf{vendredi 18 décembre 2009}.
Les rapports, au format .pdf, doivent être envoyés par email à l'adresse \textbf{galerne@cmla.ens-cachan.fr}.
Le sujet de l'email doit être \textbf{Rapport du TP3 de \textit{Nom Prénom}} et le fichier .pdf doit être intitulé \textbf{\textit{nom\_prenom}\_tp3.}
\end{description}


\subsection{Téléchargement des fichiers nécessaires}

Ouvrez un terminal et tapez y ligne par ligne les instructions suivantes pour télécharger les fichiers nécessaires et lancer Matlab (on suppose ici que le dossier \verb1tp_fourier1 a été créé dans votre dossier \texttt{home} lors du premier TP) :

\begin{verbatim}
cd tp_fourier
wget mw.cmla.ens-cachan.fr/~galerne/tp3.zip
unzip tp3.zip
cd tp3
matlab
\end{verbatim}




\section{Images numériques et polynômes trigonométriques}


Dans tout le TP on considère  une image numérique $u=(u_{m,n})$,
$(m,n)\in \llbracket 0, M-1\rrbracket \times \llbracket 0, N-1\rrbracket$,
où \textbf{les entiers $M$ et $N$ sont supposés pairs}.
La transformée de Fourier discrète (TFD) de $u$ est notée $\tilde{u}=(\tilde{u}_{m,n})$,
avec $(m,n)\in \left\llbracket -\frac{M}{2}, \frac{M}{2}-1\right\rrbracket \times
\left\llbracket -\frac{N}{2}, \frac{N}{2}-1\right\rrbracket$.
On rappelle que le polynôme trigonométrique $(M,N)$-périodique $P_u$ associé à $u$ est défini par
\begin{equation}
P_u(x,y) =
\sum_{m=-\frac{M}{2}}^{\frac{M}{2}-1} \sum_{n=-\frac{N}{2}}^{\frac{N}{2}-1}
\tilde{u}_{m,n} e^{\frac{2i\pi m x}{M}} e^{\frac{2i\pi n y}{N}}.
\label{eq_definition_Pu}
\end{equation}

Le but des TP3 et TP4 est de résoudre de manière exacte des équations aux dérivées partielles (EDP)
dans l'espace des polynômes trigonométriques.
D'un point de vue pratique il est important de se restreindre à des polynômes trigonométriques à valeur réelle.
La proposition suivante explicite une condition nécessaire et suffisante sur $u$
pour que le polynôme trigonométrique $P_u$ soit à valeur réelle
(cette proposition est l'équivalent en deux dimension de la proposition 6.4 p. 87 du polycopié).

\begin{proposition}(CNS pour que $P_u$ soit réel)\\
Le polynôme trigonométrique $P_u$ est à valeur réelle sur $\R^2$ si et seulement si
$$
\forall m\in\left\llbracket -\frac{M}{2}, \frac{M}{2}-1\right\rrbracket,~\tilde{u}_{m,-\frac{N}{2}}=0 \quad \text{et} \quad
\forall n\in\left\llbracket -\frac{N}{2}, \frac{N}{2}-1\right\rrbracket,~\tilde{u}_{-\frac{M}{2},n}=0.\\
$$
\end{proposition}

\begin{question}
Ecrivez une fonction Matlab intitulée \verb1function v = annuler_hf(u)1 qui prend en entrée une image $u$ de taille $M\times N$ et qui retourne l'image $v$ obtenue en annulant les $M+N-1$ coefficients hautes fréquences de la TFD de $u$,
de telle sorte que le polynôme triogonométrique associé à $v$ soit à valeur réelle.
Testez cette fonction sur l'image \texttt{lena.png} (ou une autre image en niveau de gris) et vérifiez que la différence entre $u$ et $v$ n'est pas ou peu perceptible.
Visualisez également la différence $u-v$ et commentez l'allure de cette fonction.
\end{question}




\section{Calcul des dérivées d'une image}

Un prérequis essentiel pour la résolution numérique d'EDP est de disposer d'algorithmes de calcul des dérivées des images.
Ici nous allons identifier l'image $u$ avec son polynôme trigonométrique associé $P_u$.
Nous allons donc implémenter des algortihmes qui calculent exactement les dérivées de $P_u$ aux points $(m,n)$.


\subsection{Gradient d'une image}

\begin{question}
En dérivant l'expression~\eqref{eq_definition_Pu} exprimez les coefficients des TFD des images
$\dfrac{\partial P_u}{\partial x}(m,n)$ et $\dfrac{\partial P_u}{\partial y}(m,n)$, $(m,n)\in \llbracket 0, M-1\rrbracket \times \llbracket 0, N-1\rrbracket$, en fonction de la TFD $\tilde{u}$ de $u$.\\
\end{question}

\begin{question}
Ecrivez une fonction Matlab
\verb1function [dxu, dyu]=gradient_tfd2(u)1
qui calcule les images $\dfrac{\partial P_u}{\partial x}(m,n)$ et $\dfrac{\partial P_u}{\partial y}(m,n)$.
Testez cette fonction sur une image sans hautes fréquences et visualisez l'image originale, les deux dérivées partielles,
ainsi que la norme du gradient en chaque point.\\
\end{question}

Ainsi calculé, le gradient est sujet aux effets de Gibbs, notament au bord de l'image.
Pour éviter ce problème on utilise une symétrisation préalable de l'image, en la rendant paire selon chaque variable (symétrisation miroir).
Cette symétrie est réalisée par la fonction \texttt{miroir} (voir le fichier \texttt{miroir.m}).\\


\begin{question}
A l'aide de votre fonction \verb1gradient_tfd21 calculez le gradient de l'image symétrisé \texttt{v = mirroir(u)}.
Extrayez de ce gradient le sous-domaine correspondant à l'image originale.
Comparez numériquement ce ``gradient symétrique'' avec le ``gradient périodique'' calculé à la question précédente.
\end{question}


\subsection{Laplacien d'une image}

On s'intéresse maintenant au calcul du laplacien d'une image en suivant la même démarche qu'à la section précédente.\\

\begin{question}
Donnez la relation entre les coefficients de la TFD de $u$ et ceux de la TFD de
$$
\Delta P_u (m,n) = \frac{\partial^2 P_u}{\partial x^2}(m,n) + \frac{\partial^2 P_u}{\partial y^2}(m,n),
~(m,n)\in \llbracket 0, M-1\rrbracket \times \llbracket 0, N-1\rrbracket.
$$
\end{question}

Comme pour le gradient nous allons implémenter deux algorithmes pour le calcul du laplacien.
Le premier va calculer le ``laplacien périodique'' $\Delta_{\operatorname{per}}u = \Delta P_u (m,n)$
qui est égale au laplacien de $P_u$.
Le deuxième va calculer le ``laplacien symétrique'' $\Delta_{\operatorname{sym}}u$, qui est égale à la restriction sur le domaine
$\llbracket 0, M-1\rrbracket \times \llbracket 0, N-1\rrbracket$ du laplacien périodique de l'image symétrisée miroir de $u$.
Plus précisément, si $\mathcal{M}(u)$ est l'image symétrisée de taille $2M\times 2N$ alors
$$
\Delta_{\operatorname{sym}}u (m,n) = \Delta_{\operatorname{per}} \left(\mathcal{M}(u)\right)(m,n),
~(m,n)\in \llbracket 0, M-1\rrbracket \times \llbracket 0, N-1\rrbracket.\\
$$


\begin{question}
Ecrivez une fonction \verb1v=laplacien_per_dft2(u)1 qui calcule le laplacien périodique de $u$.
Ecrivez ensuite une fonction \verb1v=laplacien_sym_dft2(u)1 qui calcule le laplacien symétrique de $u$.
Visualisez sur un exemple l'image originale, son laplacien périodique, son laplacien symétrique ainsi que la différence des deux laplaciens.\\
\end{question}


Remarquez que l'opération $u\mapsto \Delta_{\operatorname{per}}u$ est presque inversible puisque,
à l'exception du coefficient $\tilde{u}_{0,0}$,
chaque coefficient de Fourier est multiplié par un facteur non nul.
Ainsi seul le coefficient de Fourier $\tilde{u}_{0,0}$ n'est pas déterminé par la donné de $\Delta_{\operatorname{per}}u$
(on rappelle que $\tilde{u}_{0,0}$ représente soit la moyenne de l'image $u$ soit la somme de tous ses coefficients,
selon la convention TFD/TFD inverse).

Remarquez également que les opérateurs différentiels dépendent des conditions de bord.
Ceci a son importance en pratique.
Par exemple, pour résoudre l'équation de Poisson avec des conditions périodiques on doit utiliser l'opérateur $\Delta_{\operatorname{per}}$.
En revanche, pour résoudre l'équation de Poisson avec des conditions de Neumann (bord adiabatique) on doit utiliser l'opérateur
$\Delta_{\operatorname{sym}}$.
Pour les conditions de Dirichlet (bord isotherme), il faut utiliser un troisième laplacien obtenu à partir d'une symétrisation impaire de l'image.




\section{Composante périodique d'une image}

Dans cette section on définit et calcule la composante périodique d'une image $u$.
La composante périodique d'une image répond au problème suivant :
on cherche à déterminer une image $p = \operatorname{per}(u)$ qui soit à la fois proche de $u$
et qui soit ``analytiquement périodique'', c'est-à-dire avec de faibles sauts d'un bord à l'autre de l'image.

\begin{definition}(composante périodique)\\
Soit $u$ une image numérique.
La composante périodique $p = \operatorname{per}(u)$ est l'unique solution $p$ du problème de Poisson suivant :
\begin{equation}
\begin{cases}
\Delta_{\operatorname{per}} p = \Delta_{\operatorname{sym}} u,\\
\operatorname{moy}(p) = \operatorname{moy}(u),
\end{cases}
\label{eq_def_composante_periodique}
\end{equation}
où $\operatorname{moy}(u)$ désigne la moyenne de $u$.
\end{definition}


Le calcul de la composante périodique de $u$ s'effectue donc de la manière suivante.
\begin{enumerate}
\item Calculer $\Delta_{\operatorname{sym}} u$.
\item Calculer la TFD de $\Delta_{\operatorname{sym}} u$.
\item Diviser chaque coefficients de Fourier par le facteur multipicatif qui intervient dans le calcul de $\Delta_{\operatorname{per}} u$, à l'exception du coefficient d'indice $(0,0)$.
\item En $(0,0)$ imposer la valeur $\tilde{u}_{(0,0)}$ (sans toutefois calculer la TFD de $u$!) pour que la TFD inverse ait la même moyenne que $u$.
\item Effectuer une TFD inverse.\\
\end{enumerate}



\begin{question}
Ecrivez une fonction \verb1function p=composante_periodique(u)1
qui calcule la composante périodique d'une image (sans hautes fréquences) $u$.\\
\end{question}

\begin{question}
Calculez la composante périodique $p$ de l'image $u=$\verb1nuages_ng.png1, ainsi que la différence $s=p-u$.
Sur une même figure Matlab, visualisez les trois images $u$, $p$ et $s$ ainsi que leurs modules respectifs
(utilisez les commandes \texttt{subplot(2,3,k)} et la fonction \verb1vis_mod_img(v)1).
Commentez.\\
\label{quest_visualisation_composante_periodique}
\end{question}


% Il est possible d'analyser certains types d'image à partir du module de leur DFT, notament les images de textures.
% Par exemple, si une droite est visible dans le module de Fourier d'une texture,
% c'est que la direction orthogonale à cette droite est une direction caractéristique de la texture.\\
% 
% \begin{question}
% Reprendre la Question~\ref{quest_visualisation_composante_periodique} avec
% l'image \verb1???1 à la place de l'image \verb1nuages_ng.png1.
% En s'appuyant sur ces deux exemples,
% justifiez qu'il est plus pertinent d'utiliser la composante périodique plutôt que
% l'image originale pour effectuer une analyse spectrale de la texture.
% \end{question}



\section{Images couleur}

Dans la suite du TP nous allons utiliser des images numériques en couleur.
Une image couleur $u$ est représentée par trois images en niveau de gris $u = (u_r,u_g,u_b)$.
Chacune de ces trois images représente la quantité de rouge, de vert et de bleu en chaque pixel.
On parle alors de représentation RGB (pour red, green, blue).

Pour Matlab, une image couleur est donc un tableau de chiffre à trois dimension de taille $M\times N\times 3$.
Pour visualiser une image couleur on utilisera les fonctions \verb1vis_imgc_mat1 et \verb1vis_imgc_ext1.

A titre d'illustration, exécutez le script intitulé \verb1canaux_couleurs.m1 qui lit l'image couleur \verb1singe.png1
et affiche sur une même figure l'image couleur ainsi que ses trois canaux couleurs.


\section{Synthèse d'images par phase aléatoire}

\subsection{Bruit à phase aléatoire pour la synthèse de textures}

On a observé lors du TP1 que l'information des contours géométriques d'une image est principalement contenue dans la phase d'une image.
Inversement, l'information de la texture d'une image est plutôt contenu dans le module de Fourier d'une image.

Cette dernière observation peut être mise à profit pour synthétiser des textures.
Le bruit à phase aléatoire (BPA) associé à une image $u$ est une image aléatoire obtenu en remplaçant la phase de la TFD de $u$ par une phase aléatoire.
Comme la phase est rendue incohérente, le BPA transforme n'importe qu'elle image $u$ en une texture.
L'intérêt de ce procédé est qu'il est capable de resynthétiser certaines textures naturelles.
\\

\begin{question}
Calculez à l'aide de la fonction \verb1bpa_couleur.m1 une réalisation des BPA associés aux images
\verb1singe.png1,
\verb1mur.png1,
\verb1tissu.png1,
\verb1nuages.png1,
et \verb1briques.png1.
Parmi les quatre images de texture quelles sont celles qui semblent être bien reproduites ?\\
\label{ques_bpa}
\end{question}

Pour la question suivante nous avons besoin d'étendre la notion de composante périodique aux images couleurs.
La composante périodique d'une image couleur $u = (u_r, u_g, u_b)$ est naturellement définie comme étant l'image couleur
$p$ dont les canaux sont les composantes périodiques des canaux de $u$, c'est-à-dire
$p = (\operatorname{per}(u_r), \operatorname{per}(u_g), \operatorname{per}(u_b))$.\\

\begin{question}
Comparez le BPA associé à \verb1nuages.png1 et le BPA associé à la composante périodique de \verb1nuages.png1.
Comment expliquez-vous les différences observées ?
\end{question}



\subsection{``Fonctions typiques'' d'un espace de Sobolev}

Dans cette partie on synthétise des images qui représentent des ``fonctions typiques'' d'espace de Sobolev.
On rappelle d'abord le résultat suivant (Proposition 7.6 p. 113 du polycopié).

\begin{proposition}(norme sur les esapces $H^k_{\operatorname{per}}$)\\
Les espaces $H^k_{\operatorname{per}}\left([0,2\pi]^2\right)$ sont des espaces de Hilbert et ont pour norme équivalente
$$
\left\|f\right\|^2_{H^k_{\operatorname{per}}} = \sum_{(m,n)} |c_{(m,n)}(f)|^2 (1+m^2+n^2)^k.
$$
\end{proposition}

Les coefficients de Fourier des fonctions $f\in H^k_{\operatorname{per}}\left([0,2\pi]^2\right)$ décroissent donc rapidement.
Dans ces conditions, les coefficients de la TFD des images discrètes
$f_N : (m,n)\mapsto f\left( \frac{2\pi m}{N},\frac{2\pi m}{N}\right)$ sont de bonnes approximations des coefficients de Fourier de $f$ (en vertu du théorème 6.1 p. 89).\\

\begin{question}
Ecrivez une fonction \verb1function u = sobolev(alpha,N)1 qui calcule une image $u$ de taille $N\times N$ dont les modules des coefficients de Fourier sont
$$
\left|\tilde{u}_{m,n}\right| = \left(1+m^2+n^2\right)^{-\frac{\alpha}{2}}
$$
et dont la phase est aléatoire (utilisez la fonction \verb1phase_aleatoire1).
Choisissez trois valeurs de $\alpha$ pour visualiser une fonction $f_1$ en dehors de $H^1_{\operatorname{per}}\left([0,2\pi]^2\right)$, une fonction $f_2$ dans
$H^1_{\operatorname{per}}\left([0,2\pi]^2\right) \setminus H^2_{\operatorname{per}}\left([0,2\pi]^2\right)$,
et une fonction $f_3$ de $H^2_{\operatorname{per}}\left([0,2\pi]^2\right)$.\\
\end{question}

Comme la phase est aléatoire et que la décroissance des coefficients est choisie de manière optimale,
on peut argumenter que ces fonctions représentent la ``régularité typique'' de leurs espaces de Sobolev respectifs.

Enfin du point de vue de la synthèse d'images, remarquez que la fonction \texttt{sobolev} peut-être considérée comme un synthétiseur de nuages plus ou moins granuleux selon le paramètre $\alpha$.

% Réponse : si \alpha > k+\frac{1}{2} alors $u$ est dans $H^k_{\operatorname{per}}$
% \alpha_1 = 1 ; \alpha_2 = 1,6 ; \alpha_3 = 2,6.




\end{document}









