%matplotlib inline
from pymatbridge import Matlab
mlab = Matlab()
mlab.start()
%load_ext pymatbridge

"""
Explanation: <div align="center">
  <h2> Méthodes quantitatives en neurosciences </h2>
</div>

<div align="center">
<b><i> Cours NSC-2006, année 2015</i></b><br>
<b>Laboratoire d'introduction au traitement des signaux - Réponses</b><br>
     *Pierre Bellec, Yassine Ben Haj Ali*
</div>

Objectif:
Ce laboratoire a pour but de vous familiariser avec la notion de convolution et de système linéaire invariant dans le temps. Nous allons implémenter des convolutions avec différents noyaux, et tester la linéarité du système. Nous allons reproduire la plupart des résultats présentés dans le cours avec matlab pour des signaux temporels, et également appliquer les mêmes techniques à des images. Pour réaliser l'exercice 2 de ce laboratoire, il est nécessaire de récupérer l'<a href="http://www.imageprocessingplace.com/downloads_V3/root_downloads/image_databases/standard_test_images.zip">archive</a> suivante, et copier l'image <code>lena_gray_512.tif</code> dans votre répertoire de travail.
Ne pas tenir compte et ne pas executer cette partie du code:
End of explanation
"""


%%matlab
time_samples = 0.1*(0:1000);
bloc = zeros(size(time_samples));
bloc(300:500) = 1;
bloc = bloc + sqrt(0.01)*randn(size(time_samples));

"""
Explanation: Section 1 : Convolution de signaux temporels
<ol start="1">
  <h4><li> Commençons par générer un signal de type `bloc`, auquel on rajoute un bruit gaussien:</li></h4>
</ol>
End of explanation
"""


%%matlab
size(bloc)

"""
Explanation: Quel est la taille du signal? 

Le signal bloc est un vecteur ligne de longueur 1001:
End of explanation
"""


%%matlab
figure
plot(time_samples,bloc)
xlabel('time')
ylabel('a.u.')

"""
Explanation: La taille du bloc? 
A l’aide de la ligne bloc(300:500) = 1; , on voit que le bloc contient 201 valeurs non-nulles, avant addition
du bruit.


Quelle est la fréquence d’échantillonnage, en supposant que les
    temps d’acquisition des échantillons sont représentés par
    time_samples en secondes? 
A l’aide de la ligne
time_samples = 0.1*(0:1000); on voit que les échantillons
temporels sont séparés de 0.1 s, soit une fréquence
d’échantillonnage de 1 Hz.


Quelle est la durée d’acquisition?
L’enregistrement dure 100
secondes.


Générez un graphe représentant bloc avec l’axe temporel
    adéquat (voir la fonction plot).
Pour générer un graphe avec
l’axe temporel adéquat, on utilise plot avec deux arguments, ainsi que les commandes xlabel et ylabel     pour labéliser les axes.
End of explanation
"""


%%matlab
noyau = ones(1,3);
noyau = noyau/sum(abs(noyau));
bloc_conv = conv(bloc,noyau,'same');

"""
Explanation: <ol start="2">
  <h4><li> On va maintenant construire un noyau, et appliquer la convolution sur le signal `bloc`</li></h4>
</ol>
End of explanation
"""


%%matlab
figure
stem(noyau)
axis([0 4 0 0.5])

"""
Explanation: Faire un graphe avec le noyau (note: stem est plus adapté que plot pour les signaux courts, et \mcode{axis} peut être utilisé pour ajuster les limites des axes).
End of explanation
"""


%%matlab
figure
plot(time_samples,bloc,'b')
hold on
plot(time_samples,bloc_conv,'r')

"""
Explanation: Quelle opération est associée avec la convolution basée sur ce noyau? 
Ce noyau implémente une moyenne mobile.


Superposez les signaux bloc et bloc_conv sur un même graphe. Qu'observez-vous?
End of explanation
"""


%%matlab
for nn = [7 100 300]
    figure
    noyau = ones(1,nn); noyau = noyau/sum(abs(noyau)); bloc_conv = conv(bloc,noyau,'same');
    hold on, plot(time_samples,bloc,'b'), plot(time_samples,bloc_conv,'r')
end

"""
Explanation: On observe une diminution de la variance du bruit, lié au moyennage.

<ol start="3">
  <h4><li> Répétez la question précédente en utilisant un noyau de longueur 7, puis 100, puis 300.</li></h4>
</ol>

Qu'observez-vous?
End of explanation
"""


%%matlab
noyau = [-ones(1,10) 0 ones(1,10)];
noyau = noyau/sum(abs(noyau));
bloc_conv = conv(bloc,noyau,'same');

"""
Explanation: On observe que pour des filtres de taille croissante, le bruit devient très faible, mais que les transitions du bloc deviennent de moins en moins abruptes. Pour un filtre dépassant la taille du bloc (300), la valeur maximale du bloc aprés convolution devient plus petite que la valeur moyenne dans le bloc avant convolution.

<ol start="4">
  <h4><li> On va maintenant s'intéresser à un nouveau noyau, que l'on va appliquer sur le même signal `bloc`.</li></h4>
</ol>
End of explanation
"""


%%matlab
figure
stem(noyau)
axis([0 22 -0.06 0.06 ])

"""
Explanation: Faire un graphe avec le noyau.
End of explanation
"""


%%matlab
figure
plot(time_samples,bloc,'b')
hold on
plot(time_samples,bloc_conv,'r')

"""
Explanation: Quelle opération est associée avec la convolution basée sur ce noyau? 
Ce noyau est sensible aux variations entre échantillons temporels voisins.


Superposez les signaux bloc et bloc_conv sur un même graphe. Qu'observez-vous?
End of explanation
"""


%%matlab
time_samples = 1*(0:100);
bloc = zeros(size(time_samples));
bloc(30:38) = 1;

"""
Explanation: Le signal convolué comporte deux pics importants, qui coincident avec les bords du bloc. une transition de la ligne de base vers le bloc induit une réponse négative, tandis qu'une transition du bloc vers la ligne de base génère une réponse positive.

<ol start="5">
  <h4><li>  On va maintenant construire un signal bloc plus court, et sans bruit, pour illustrer le concept de réponse à une impulsion finie. </li></h4>
</ol>
End of explanation
"""


%%matlab
temps_pic = 3;
noyau = [linspace(0,1,temps_pic+1) linspace(1,-0.3,temps_pic/2) linspace(-0.3,0,temps_pic/2)];
noyau = [zeros(1,length(noyau)-1) noyau];

"""
Explanation: On va également construire un noyau qui ressemble à une réponse hémodynamique, simplifiée:
End of explanation
"""


%%matlab
figure
stem(noyau)

"""
Explanation: Faire un graphe avec le noyau.
End of explanation
"""


%%matlab
bloc_conv = conv(bloc,noyau,'same');

"""
Explanation: Quelle opération est associée avec la convolution avec le noyau? 
Ce noyau est un modèle grossier d'une fonction de réponse hémodynamique cérébrale.


<ol start="6">
  <h4><li>  On va maintenant appliquer la convolution avec le noyau au signal `bloc`: </li></h4>
</ol>
End of explanation
"""


%%matlab
figure
plot(time_samples,bloc,'b-o')
hold on
plot(time_samples,bloc_conv,'r-o')

"""
Explanation: Superposez bloc et bloc_conv sur un même graphe.
End of explanation
"""


%%matlab
bloc2 = 2*bloc;

"""
Explanation: Notez l'option d'affichage pour ajouter un marqueur à chaque point de mesure.



Combien de temps met la réponse à atteindre son pic? 
On voit que trois points séparent le début du bloc du pic de la réponse. En regardant la définition de time samples, on voit que la fréquence d'échantillonnage est de 1 Hz. La réponse et donc 3 secondes à atteindre son pic.


<ol start="7">
  <h4><li> Construire un signal similaire à bloc, mais avec une amplitude du bloc deux fois plus grande. </li></h4>
</ol>
End of explanation
"""


%%matlab
bloc_conv2 = conv(bloc2,noyau,'same');
plot(time_samples,bloc_conv-bloc_conv2,'b-o')

"""
Explanation: <ol start="8">
  <h4><li> Comparez la réponse au gros bloc, avec le double de la réponse au petit bloc, en faisant un graphe de la différence des deux signaux.</li></h4>
</ol>
End of explanation
"""


%%matlab
impulsion = zeros(size(bloc));
impulsion(50) = 1;
impulsion_conv = conv(impulsion,noyau,'same');
figure
plot(time_samples,impulsion,'b-o')
hold on
plot(time_samples,impulsion_conv,'r-o')

"""
Explanation: <ol start="9">
  <h4><li> Est ce que la propriété de mise à l'échelle est respectée?</li></h4>
</ol>

On voit que la différence entre `bloc` et `bloc2` est nulle, les deux signaux sont donc identiques. La propriété de mise à l'échelle semble être ici respectée.

<ol start="10">
  <h4><li> Construisez un nouveau stimulus `impulsion`, contenant une seul impulsion unitaire à l'échantillon numéro 50. Calculez la réponse du système à ce nouveau stimulus, appelée `impulsion_conv`, et représentez-la sur un graphe. </li></h4>
</ol>
End of explanation
"""


%%matlab
somme = bloc+impulsion;
reponse_somme = conv(somme,noyau,'same');
figure
plot(time_samples,somme,'b-o')
hold on
plot(time_samples,reponse_somme,'r-o')

"""
Explanation: <ol start="11">
  <h4><li> Calculer la réponse à la somme de `bloc` et `impulsion`, et représentez la sur un graphe.</li></h4>
</ol>
End of explanation
"""


%%matlab
somme_reponse = impulsion_conv + bloc_conv;
figure
plot(time_samples,reponse_somme - somme_reponse,'b-o')

"""
Explanation: <ol start="12">
  <h4><li> Comparez la réponse à `bloc+impulsion` avec la somme de `bloc_conv` et `impulsion_conv`, en faisant un graphe de la différence des deux signaux. La propriété d'additivité est-elle respectée?</li></h4>
</ol>
End of explanation
"""


%%matlab
impulsion2 = zeros(size(bloc));
impulsion2(70) = 1;
impulsion2_conv = conv(impulsion2,noyau,'same');
figure
plot(time_samples,impulsion_conv,'r-o')
hold on
plot(time_samples,impulsion2_conv,'b-o')

"""
Explanation: La réponse à la somme des signaux est identique à la somme des réponses. La propriété d'additivité est ici respectée.

<ol start="13">
  <h4><li> Répétez la question 8 avec une impulsion localisé à l'échantillon 70. Est ce que la propriété d'invariance dans le temps est respectée?</li></h4>
</ol>
End of explanation
"""


%%matlab
img = imread ('lena_gray_512.tif');
img = double(img);

"""
Explanation: On observe que la réponse à `impulsion2` est identique à `impulsion`, décalée dans le temps. La propriété d'invariance temporelle semble respectée.

Section 2 : Convolution d'images

En premier lieu, il faut charger l'image lena_gray_512.tif à l'aide de la fonction imread():
End of explanation
"""


%%matlab
figure
imagesc(img,[0 255]) % le 2ieme argument ajuste le minimum/maximum 
colorbar % Ajoutez une echelle de couleur
colormap gray % Utilisez une echelle en niveau de gris

"""
Explanation: La fonction imagesc() nous permet de visualiser une image binaire comme suit:
End of explanation
"""


%%matlab
whos img

"""
Explanation: Quelle est la taille de cette image? son type? Combien y-a-t-il de niveaux de gris dans cette image?
End of explanation
"""


%%matlab
noyau1 = ones(5,5);
noyau1 = noyau1/sum(sum(abs(noyau1)));

"""
Explanation: l'image et de dimention 250x250 et de type binaire, elle a 250 niveau de gris


Nous allons maintenant créer un noyau pour la convolution, qui est une autre petite image définissant les poids d'une moyenne pondérée dans chaque petit carré (ici, 5x5) de l'image:
End of explanation
"""


%%matlab
imagesc(noyau1)
colorbar
colormap gray

"""
Explanation: Visualisez le noyau avec la commande imagesc.
End of explanation
"""


%%matlab
img1 = conv2(img,noyau1,'same');

"""
Explanation: Appliquons la convolution entre img et le noyau:
End of explanation
"""


%%matlab
figure
imagesc(img,[0 255]); colormap gray; colorbar
figure
imagesc(img1,[0 255]); colormap gray; colorbar

"""
Explanation: Visualisez les deux images, avant et après la convolution, en prenant soin d'utiliser les mêmes valeurs min et max dans les deux cas. Quelle opération implémente ce noyau? Que remarquez vous sur l'image convoluée?
End of explanation
"""


%%matlab
noyau2 = ones(20,20);
noyau2 = noyau2/sum(sum(abs(noyau2)));
img2 = conv2(img,noyau2,'same');
figure
imagesc(img,[0 255]); colormap gray; colorbar
figure
imagesc(img2,[0 255]); colormap gray; colorbar

"""
Explanation: Ce noyau implement un moyenne mobile, le resultat sur l'image est une operation de lissage.
End of explanation
"""


%%matlab
noyau3 = [ -ones(5,2) zeros(5,1) ones(5,2) ]
imagesc (noyau3); colorbar; colormap gray

"""
Explanation: Le lissage est plus agressif(floutage), il y a perte de cerataines transitions de contrast dans l'image.


Construisons maintenant un nouveau noyau:
End of explanation
"""


%%matlab
img3 = conv2(img,noyau3,'same');
figure
imagesc(img,[0 250]); colormap jet; colorbar
figure
imagesc(img3,[0 250]); colormap jet; colorbar

"""
Explanation: Appliquez la convolution sur l'image img avec le nouveau noyau noyau3 et enregistrez la nouvelle image convoluée dans la variable img3. Visualisez les images avant/après convolution, en utilisant l'échelle de couleur jet. Que remarquez vous?
End of explanation
"""


%%matlab
noyau4 = noyau3'
figure
imagesc (noyau4); colorbar; colormap gray

%%matlab
img4 = conv2(img,noyau4,'same');
figure
imagesc(img,[0 250]); colormap jet; colorbar
figure
imagesc(img4,[0 250]); colormap jet; colorbar

"""
Explanation: Ce noyau permet une detection des bords veticaux dans l'image


Répétez les deux questions précédentes avec ces autres noyaux:

1- noyau4:
End of explanation
"""


%%matlab
noyau5 = tril(-ones(size(noyau4))) + triu(ones(size(noyau4)))
figure
imagesc (noyau5)
colorbar
colormap gray

%%matlab
img5 = conv2(img,noyau5,'same');
figure
imagesc(img,[0 250]); colormap jet; colorbar
figure
imagesc(img5,[0 250]); colormap jet; colorbar

"""
Explanation: C'est un noyau qui detecte le bords horizontaux de l'image
2- noyau5:
End of explanation
"""


%%matlab
img_sum = abs(img3)+abs(img4)+abs(img5);
figure
imagesc(img_sum,[0 250]); colormap jet; colorbar

"""
Explanation: Ce noyau detecte les bords en ongle de ~45 degré


Soit img4 (respectivement img5) le résultat de la convolution avec noyau4 (respectivement noyau5). Visualisez et interprétez l'image abs(img3)+abs(img4)+abs(img5).
End of explanation
"""


