UP | HOME

Introduction à la programmation Octave/Matlab

Table des matières

Responsable

Arnaud Malapert.

Introduction

Ce cours s'étend sur 6 séances de l'U.E "Informatique" du Master Pro ISAB (Imagerie et Systèmes Appliqués en Biologie). La partie pratique est une initiation au logiciel Gnu Octave. Ce logiciel est un clone libre du logiciel Matlab très répandu dans l'industrie. Gnu Octave permet d'effectuer des calculs (matriciels, statistiques, analytiques, etc), de la visualisation avec Gnuplot (galerie) et du traitement d'images.
Octave/Matlab est un langage interprété de haut niveau pour le calcul numérique. Il permet de résoudre numériquement des problèmes linéaires ou non linéaires et de réaliser toutes sortes d'expérimentation numériques, ou des analyses statistiques. Il offre aussi des fonctionnalités graphiques avancées pour la visualisation et le traitement d'images. Octave est utilisé soit en ligne de commande de manière interactive ou pas. Les langages Octave et Matlab sont très similaires et il est possible d'écrire des programmes portables.
Octave/Matlab est aussi très complémentaire d'un tableur (MS Excel, OO Calc, etc):

  • Octave : automatisation de calculs complexes ; tracer de meilleurs graphiques.
  • Tableur : lecture et formatage des tableaux de résultats.

Contenu

  • Les transparents du cours sont disponibles ici.
  • Les notes de cours sont disponibles ici.

Le contenu des transparents et des notes de cours est identique.

Ressources électroniques

À propos de la machine virtuelle Linux

  1. Installer le logiciel VirtualBox.
  2. Importer la machine virtuelle (virtual machine - VM) fournie en cours.
  3. Démarrer la VM
  4. Bienvenue sur xubuntu, une variante d'ubuntu.
  5. Quelques logiciels utiles.
  6. J'ai le reflexe de chercher de l'aide, notamment sur les forums.
  7. Démarrer Octave :
    • Démarrer un terminal et tapez octave. Lisez la partie 2 du guide (plutôt que nano, utiliser gedit ou emacs) ;
    • ou démarrer l'environnement de développement intégré QtOctave.

Projet

Mise en application de techniques basiques de programmation en nous appuyant sur les cours 1 et 2 du cours de traitement d'images de Y. Tarabalka.

Nous supposerons que les images sont en noir et blanc et que l'intensité est comprise entre 0 et 2**bitDepth-1bitDepth est un nombre de bits.

Utiliser les fonctions imread pour importer une image et imwrite pour exporter nos résultats. L'image de test isab/aaa.png est disponible ici.

# WARNING ! The file must exist.
infile = 'isab/aaa.png';
bitDepth = imfinfo(infile).BitDepth
X = imread(infile);

Fonctions basiques (cours 1)

  1. Programmer une fonction imgNegative(X, bitDepth) renvoyant le négatif d'une image X.
  2. Programmer une fonction imgThreshold(X, threshold, bitDepth) renvoyant l'image X après l'application d'une fonction de seuil.
  3. Programmer une fonction imgQuantify(X, newBitDepth, bitDepth) avec des intensités codées sur newBitDepth bits.
  4. Programmer une fonction imgSample(X, newWidth, newHeight) redimensionnant l'image par un échantillonnage simple, on sélectionne les pixels en quadrillant l'image.
function Y = imgNegative(X, bitDepth = 8) 
  Y = 2**bitDepth - 1 - X;
endfunction

function Y = imgThreshold(X, threshold, bitDepth = 8) 
Y =  (X >= threshold) * (2**bitDepth-1);
endfunction

function Y = imgQuantify(X, newBitDepth, bitDepth) 
  Y =  X - mod( X, 2**(bitDepth-newBitDepth));
endfunction

function Y = imgSample(X, newWidth, newHeight)
  Y =  X( round(linspace(1, size(X,1), newWidth)), round(linspace(1, size(X,2), newHeight)) );
endfunction
outfile = 'isab/aaa_utils.png';
Y = [ X imgNegative(X, bitDepth) imgThreshold(X, 127, bitDepth) imgQuantify(X, 3, 8) ];
imwrite(Y,  gray(256), outfile)

aaa_utils.png

outfile = 'isab/aaa_reduced.png';
Y = imgSample(X,128,128);
imwrite(Y,  gray(256), outfile)

aaa_reduced.png

Génération de bruits aléatoires (p. 4-6, cours 2)

  1. Programmer la fonction imgPattern(n, bitDepth) renvoyant le motif de test.
  2. Programmer une fonction imgGNoise(X, mean, variance, bitDepth) renvoyant l'image X avec un bruit gaussien.
  3. Programmer une fonction imgUNoise(X, a, b, bitDepth) renvoyant l'image X avec un bruit uniforme.
  4. Programmer une fonction imgSPNoise(X, probSP, proba, a, b) renvoyant l'image X avec un bruit poivre et sel.
function Y = imgPattern(n = 300, bitDepth = 8)
  intens = idivide(2**bitDepth-1, 6);
  s = idivide(n, 6);
  Y = intens * ones(n); 
  IND = s+1:(n-s);
  Y(IND,IND) = 3*intens;
  IND = 2*s+1:(n-2*s); 
  Y(IND,IND) = 5*intens;
endfunction

function Y = imgGNoise(X, mean = 0, variance = 5, bitDepth = 8)
  Y = round(X + normrnd( mean, sqrt(variance), size(X)));
  Y = max(0, min(Y, 2**bitDepth-1));
endfunction

function Y = imgUNoise (X, a = -5, b = 5, bitDepth = 8)
 Y = X + randi([a b], size(X));
 Y = max(0, min(Y, 2**bitDepth-1));
endfunction

function Y = imgSPNoise(X, probSP = 0.05, proba = 0.5, a = 0, b =255)
  Y = X;
  IND = find(rand(size(X)) <= probSP);
  SP = rand(1, length(IND) ) <= proba;
  Y(IND(SP)) = a;
  Y(IND(!SP)) = b;
endfunction
outfile = 'isab/pattern.png';
P = imgPattern(180);
I = [P imgGNoise(P, 0, 20) imgUNoise(P, -20, 20) imgSPNoise(P, 0.1, 0.5, 0, 2**8-1) ];
imwrite(I,  gray(256), outfile)

pattern.png

Convolutions (p. 14-20, cours 2)

  1. Programmer une fonction imgZeroPadding(X, MASK) renvoyant l'image X étendu par zero-padding par application d'un masque.
  2. Programmer une fonction imgLinearFilter(X, MASK) renvoyant l'image X après lissage par un filtre linéaire (p. 14, cours 2).
  3. Programmer une fonction imgMeanFilter(X, radius) qui applique un filtre par moyenne en spécialisant imgLinearFilter.
  4. Programmer une fonction imgCrossFilter(X, radius) en spécialisant imgLinearFilter.
  5. Programmer une fonction imgMedianFilter(X, radius) renvoyant l'image X après lissage par un filtre médian.
  6. Programmer une fonction imgMiddleFilter(X, radius) renvoyant l'image X après lissage par un filtre milieu.
function Y = imgZeroPadding(X, MASK)
  Y = X;
  for i= 1:length(size(X)) 
    padlen = idivide(size(MASK, i), 2);
    Y = prepad(Y, size(Y, i) + padlen, 0, i);
    Y = postpad(Y, size(Y, i) + padlen, 0, i);
  endfor
endfunction


function Y = imgLinearFilter(X, MASK)
  Y = zeros(size(X));
  X = imgZeroPadding(X, MASK);
  di = idivide(size(MASK,1), 2);
  dj = idivide(size(MASK,2), 2);
  I = 1+di:(size(X, 1)-di);
  J = 1+dj:(size(X, 2)-dj);
  for i = I
    for j = J
      VX = X(i-di:i+di, j-dj:j+dj);
      Y(i-di, j-dj) = round(sum(sum(VX.*MASK))/sum(sum(MASK)));
    endfor
  endfor
endfunction

function Y = imgMeanFilter(X, radius = 1)
  Y = imgLinearFilter(X, ones(2*radius+1));
endfunction

function Y = imgCrossFilter(X, radius = 1)
  MASK = ones(2*radius+1);
  MASK(radius+1, :) = 2;
  MASK(:, radius+1) = 2;
  MASK(radius+1, radius+1) = 4;
  Y = imgLinearFilter(X, MASK);
endfunction


function Y = imgOrderFilter(X, MASK, useMiddle = 0)
  Y = zeros(size(X));
  X = imgZeroPadding(X, MASK);
  di = idivide(size(MASK,1), 2);
  dj = idivide(size(MASK,2), 2);
  I = 1+di:(size(X, 1)-di);
  J = 1+dj:(size(X, 2)-dj);
  for i = I
    for j = J
      VX = vec( X(i-di:i+di, j-dj:j+dj) );
      if(useMiddle == 0)
        Y(i-di, j-dj) = median(VX);
      else
        Y(i-di, j-dj) = idivide(max(VX)+min(VX), 2);
      endif
    endfor
  endfor
 endfunction

function Y = imgMedianFilter(X, radius = 1)
  Y = imgOrderFilter(X, ones(2*radius+1), 0);
endfunction

function Y = imgMiddleFilter(X, radius = 1)
  Y = imgOrderFilter(X, ones(2*radius+1), 1);
endfunction
  • Exemple de zero-padding
    imgZeroPadding(imgPattern(6), ones(3))
    
    ans =
    
         0     0     0     0     0     0     0     0
         0    42    42    42    42    42    42     0
         0    42   126   126   126   126    42     0
         0    42   126   210   210   126    42     0
         0    42   126   210   210   126    42     0
         0    42   126   126   126   126    42     0
         0    42    42    42    42    42    42     0
         0     0     0     0     0     0     0     0
    
  • Lisser une image avec les filtres moyens, en croix, médian, et milieu (colonnes) de forme carré et de rayons radius (lignes).
    function I = applyFilters(X, radius = 1, bitDepth = 8)
      I = [];
      for r=radius
        Y = postpad(X, size(X, 2) + 10, 0, 2);
        Y = [ Y imgMeanFilter(X, r)];
        Y = postpad(Y, size(Y, 2) + 10, 0, 2);
        Y = [ Y imgCrossFilter(X, r)];
        Y = postpad(Y, size(Y, 2) + 10, 0, 2);
        Y = [ Y imgMedianFilter(X, r)];
        Y = postpad(Y, size(Y, 2) + 10, 0, 2);
        Y = [ Y imgMiddleFilter(X, r)];
        Y = postpad(Y, size(Y, 1) + 10, 0, 1);
        I = [I ; Y];
      endfor
    endfunction
    
    
    P = imgPattern(180);
    cmap = gray(256);
    radius = 1:3;
    
  • Lisser le motif avec un bruit gaussien avec les filtres moyens, en croix, médian, et milieu (colonnes) de forme carré et de rayons de 1 à 3 (lignes).
    outfile = 'isab/gaussian_noise_filtering.png';
    B = imgGNoise(P, 0, 20);
    imwrite(applyFilters(B, radius),  cmap, outfile)
    

    gaussian_noise_filtering.png

  • Lisser le motif avec un bruit uniforme.
    outfile = 'isab/uniform_noise_filtering.png';
    B = imgUNoise(P, -20, 20);
    imwrite(applyFilters(B, radius),  cmap, outfile)
    

    uniform_noise_filtering.png

  • Lisser le motif avec un bruit poivre et sel.
    outfile = 'isab/sp_noise_filtering.png';
    B = imgSPNoise(P, 0.1,  0.5, 0, 2**8-1);
    imwrite(applyFilters(B, radius),  cmap, outfile)
    

    sp_noise_filtering.png

Rendu du projet

Le projet se réalise en binôme et sera rendu dans une boîte de dépôt JALON le <2015-12-16 mer.>. Le rendu est composé d'un seul fichier nom1_nom2.m contenant le code source commenté. Le fichier doit contenir toutes les fonctions du projet.

  • Évaluation du projet

    Le fichier main.m exécute une demonstration du projet (lien). Attention, il est essentiel que le fichier main.m s'exécute sans erreur. En supposant que les deux fichiers soient dans le répertoire de travail.

    source("nom1_nom2.m")
    source("main.m")
    
  • Fonction non implémentée

    Une fonction non implémentée doit renvoyer l'image X sans modification et afficher un message.

    function Y = imgNegative(X, bitDepth = 8) 
      warning("Function imgNegative is not implemented.");
      Y = X;
    endfunction
    imgNegative( zeros(10) );
    
    octave> warning: Function imgNegative is not implemented.
    

Created: 2017-08-02 mer. 12:08