Discussion :

[tex50]

Bonjour

Voila je dois réaliser un exercice sous MATLAB, mais je suis novice et j'ai besoin de votre aide.

J'ai des soucis à plusieurs questions , notamment à la 8) et 9)

J'aimerais qu'on m'aide à réaliser cet exercice.
Donc si vous avez le temps, merci de vérifier mon code et de m'aider à le compléter. C'est très gentil de votre part.


Initialisation de l'environnement

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
clear all; %  Supprime les variables de la mémoire
close all; %  Ferme l'ensemble des figures
clc; %  Efface les écrits en fenètre

1. Lire l'image du fantôme mesuré ('fantome_mesure.bmp')
   

   Code :

   Sélectionner tout - Visualiser dans une fenêtre à part

   1
   2
   3
   I_fantome=imread('fantome_mesure.bmp'); % la fonction "imread" permet de lire le fichier "fantome_mesure.bmp"
   figure(1);
   subplot(111);imshow(I_fantome);

2. Fabriquer le sinogramme R de l'image du fantôme mesuré en utilisant 180 projections de 1°
   

   Code :

   Sélectionner tout - Visualiser dans une fenêtre à part

   1
   2
   theta1=[0:1:180]
   [R,xp] = radon(I_fantome,theta1);

3. Fabriquer Rnorm le sinogramme R normalisé par rapport à sa valeur max
   

   Code :

   Sélectionner tout - Visualiser dans une fenêtre à part

   1
   2
   theta2=[0:180]
   [R_norm,xp] = radon(R,theta2);

4. Ajouter à Rnorm un bruit gaussien de moyenne 0 et de variance 0.005 pour obtenir Rb, sinogramme bruité
   

   Code :

   Sélectionner tout - Visualiser dans une fenêtre à part

   Rb = imnoise(R_norm,'gaussian',0,0.005);

5. Fabriquer un filtre gaussien de dimension 5x5 et d'écart-type 0.8 puis l'appliquer à Rb pour produire une version filtrée Rb_f
   

   Code :

   Sélectionner tout - Visualiser dans une fenêtre à part

   1
   2
   filtregauss = fspecial('gaussian',[5 5] , 0.8);
   Rb_f = imfilter(Rb,filtregauss);

6. Fabriquer les images reconstruites IR, IRb et IRb_f correspondantes par transformée de Radon inverse (interpolation linéaire, fenêtre Hamming)
   

   Code :

   Sélectionner tout - Visualiser dans une fenêtre à part

   1
   2
   3
   IR = iradon(R,0:1:180,'linear','hamming');
   IRb= iradon(Rb,0:180,'linear','hamming');
   IRb_f = iradon(Rb_f,0:180,'linear','hamming');

7. Déterminer les paramètres de transformation (translations uniquement) par recherche de la corrélation max avec une image fantôme de référence 256x256
   

   Code :

   Sélectionner tout - Visualiser dans une fenêtre à part

   1
   2
   3
   4
   5
   6
   7
   8
   9
   10
   11
   12
   13
   14
   15
   16
   17
   18
   19
   20
   21
   22
   23
   24
   25
   26
   27
   28
   29
   30
   %translation 
   X=100;
   Y=50;
   T=[1 0 0;0 1 0;X Y 1];
   tform = maketform('affine',T);
    
   phantom_translated=imtransform(I_fantome,tform,'bicubic',...
                      'XData', [1 (256+X)],... 
                      'YData', [1 (256+Y)]);
    
   Correlation_max=normxcorr2(I_fantome,phantom_translated);
    
   [tX tY] = size(Correlation_max);
    
   max=Correlation_max(1,1)
   for i=1:(tY-1)
       for j=1:(tX-1)
           if(Correlation_max(j,i)>max)
               max=Correlation_max(j,i);
               mX=j;
               mY=i;
           end
       end
   end
    
    
   figure(2)
   subplot(131);imshow(I_fantome);title('image originale');
   subplot(132);imshow(phantom_translated);title('translation');
   subplot(133);imshow(Correlation_max);title('Corrélation max');

8. Calculer l'erreur entre l'image fantôme initial et chacune des 3 images reconstruites (de même taille) comme la somme des différences au carré pixel à pixel
   
   
   
9. Calculer les rapports S/N en dB = 20log(somme(fantome)/erreur)
   
   
   
10. Afficher sur la même figure R, Rb et Rb_f et IR, IRb et IRb_f. Ajouter dans les légendes les valeurs des erreurs calculées.
    

    Code :

    Sélectionner tout - Visualiser dans une fenêtre à part

    1
    2
    3
    4
    5
    6
    7
    figure(3);
    subplot(331);imshow(R);title('Sinogramme R');
    subplot(332);imshow(Rb);title('Sinogramme bruité par un bruit gaussien');
    subplot(333);imshow(Rb_f);title('Sinogramme filtré par un bruit gaussien');
    subplot(334);imshow(IR);title('Image IR reconstruite ');
    subplot(335);imshow(IRb);title('Image IRb reconstruite ');
    subplot(336);imshow(IRb_f);title('Image IRbf reconstruite ');

[Jean Dumoncel]

Bonjour,

1. Dans :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

theta1=[0:1:180]

les crochets et le 1 sont inutiles, par défaut le pas est de 1 :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

theta1=0:180;

Petite remarque : on te demande d'utiliser 180 projections, le vecteur serait plutôt 0:179 qui est la valeur par défaut de la fonction radon.

3. tu fais exactement le même calcul que pour la question 2 : tu n'as pas normalisé en divisant par le maximum de R_norm. Attention au format des données de ton image : si I_fantome est au format uint8, il faudra d'abord la convertir en double avant d'effectuer la division.

7. pour la recherche du maximum, premièrement, il faut éviter d'utiliser un nom de variable qui correspond à une fonction matlab : max est déjà une fonction matlab. Ensuite :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
max=Correlation_max(1,1)
for i=1:(tY-1)
    for j=1:(tX-1)
        if(Correlation_max(j,i)>max)
            max=Correlation_max(j,i);
            mX=j;
            mY=i;
        end
    end
end

peut s'écrire justement en utilisant la fonction max :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
[maximum mI] = max(Correlation_max(:));
[mX mY] = ind2sub(size(Correlation_max),mI);

8. Pour calculer l'erreur, ce n'est pas très compliqué, il suffit de faire la différence entre tes 2 images, puis mettre chaque élément au carré (à l'aide de l'opérateur ".^2") puis de faire la somme).

[tex50]

Merci beaucoup Magelan pour cette réponse.

Pour la question 3: mon image est au format .bmp
J'ai fais comme dans la question 2 car je n'arrive pas à normaliser. Peux tu me donner la démarche à suivre.

Désolé, je suis novice sur Matlab, je commence le traitement d'images cette année.

Pourrez tu me décrire d'avantage pour la 8) et la 9).

Merci beaucoup de prendre le temps de me répondre.

[Jean Dumoncel]

Pour la 3, il suffit de faire :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

R_norm = R/max(R(:));

Oublie mon histoire de format, car on normalise la transformée de radon et non pas l'image elle-même.

Pour la 8, je ne vois pas quoi rajouter :
il suffit de faire la différence entre tes 2 images, puis mettre chaque élément au carré (à l'aide de l'opérateur ".^2") puis de faire la somme).
La différence entre les 2 images au moins, tu arrives à la faire?

[tex50]

Merci encore pour ta réponse.

Pour la 8) je dois faire comme ceci?

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
Erreur1 = I_fantome - phantom_translated
Erreur2 = Erreur1^2
Erreur3 = I_fantome - Correlation_max
Erreur4 = Erreur3^2
Erreur5 = Erreur2 + Erreur4

J'ai un message d'erreur car les matrices doivent être de la même taille pour les opérations.

Peux-tu encore m'aider s'il te plaît.
Désolé pour mon très faible niveau

[Invité]

Bonjour,

Je pense pourtant que magelan t'a plutôt bien aidé et a même devancé cette question

il suffit de faire la différence entre tes 2 images, puis mettre chaque élément au carré (à l'aide de l'opérateur ".^2") puis de faire la somme).
La différence entre les 2 images au moins, tu arrives à la faire?

Dès que l'on obtient une erreur que l'on ne comprends pas, premier réflexe: Messages d'erreur et avertissements
Ainsi donc Matrix must be square.

Ensuite pour pouvoir faire la soustraction entre deux images, il faudrait en effet avoir deux images de même taille
Dans ce sens je te conseille de revoir ta ligne

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
phantom_translated=imtransform(I_fantome,tform,'bicubic',...
                   'XData', [1 (256+X)],... 
                   'YData', [1 (256+Y)]);

[tex50]

Merci pour vos messages,

J'ai beau revoir toutes les lignes, je ne vois où est l'erreur.
Je ne vois donc pas la procédure à suivre pour soustraire les images.

Est-ce que quelqu'un pourrez m'écrire le code pour la question 8??

Merci d'avance.

[Jean Dumoncel]

Calculer l'erreur entre l'image fantôme initial et chacune des 3 images reconstruites (de même taille) comme la somme des différences au carré pixel à pixel

Ensuite pour pouvoir faire la soustraction entre deux images, il faudrait en effet avoir deux images de même taille
Dans ce sens je te conseille de revoir ta ligne

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
phantom_translated=imtransform(I_fantome,tform,'bicubic',...
                   'XData', [1 (256+X)],... 
                   'YData', [1 (256+Y)]);

Quand tu fais la translation, imtransform génère une image de taille supérieur à I_fantome, tu peux le voir en faisant :
size(phantom_translated)
Donc il faut couper phantom_translated, par exemple en faisant :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

phantom_translated = phantom_translated(1:256,1:256);

[Invité]

Et si X et/ou Y est négatif?
Pour cette raison que je me concentrerais plutôt sur la ligne contenant IMTRANSFORM.