R31 – Initiation au

traitement mathématique d’images

avec Matlab/Octave
L2 MATH & PC & SI - parcours renforcé

Recueil d’exercices corrigés et aide-mémoire.

Licence Sciences et Techniques

Gloria Faccanoni
i https://moodle.univ-tln.fr/course/view.php?id=5361

Année 2023 – 2024

Dernière mise-à-jour : Lundi 28 août 2023

Table des matières

1 Introduction à Octave/Matlab 3
1.1 Les environnements MATLAB et Octave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Installation(s) et version(s) en ligne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Premiers pas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Notions de base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.5 Commentaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.6 Affichage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.7 Opérations arithmétiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.8 Division euclidienne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.9 Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.10 Fonctions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.11 Graphes de fonctions R → R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.12 Structures itératives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.13 Structure conditionnelle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.14 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2 Traitement mathématique des images numériques 31

2.1 Introduction : lecture, affichage, sauvegarde . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.2 Manipulations élémentaires : déplacements des pixels . . . . . . . . . . . . . . . . . . . . . . 39
2.3 Manipulations élémentaires : modification des valeurs des pixels . . . . . . . . . . . . . . . . 46
2.4 Détection des bords . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
2.5 Masques (filtrage par convolution) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
2.6 Images à couleurs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
2.7 Décomposition en valeurs singulières et compression JPG . . . . . . . . . . . . . . . . . . 66
Ce fascicule est un support pour le cours d’initiation au traitement numérique d’images avec Matlab/Octave de la deuxième
année d’une Licence Scientifique.

R32 Modélisation physique – Traitement numérique d’images

CM-TP 12h 4 séances de 3h

Séance 1 : TP : Rappels Matlab/Octave

Séance 2 : TP : concepts de base, manipulations élémentaires, résolution, [transformation du photomaton et du boulanger]
Séance 3 : TP : statistiques, modification du contraste, quantification
Séance 4 : TP : débruitage, détection des bords, images à couleur
CC

Gloria FACCANONI

IMATH Bâtiment M-117 T 0033 (0)4 83 16 66 72

Université de Toulon
Avenue de l’université B [email protected]
83957 LA GARDE - FRANCE i http://faccanoni.univ-tln.fr
Chapitre 1
Introduction à Octave/Matlab
Nous illustrerons les concepts vu en cours à l’aide de MATLAB (MATrix LABoratory), un environnement de programmation et
de visualisation. Nous utiliserons aussi GNU Octave (en abrégé Octave) qui est un logiciel libre distribué sous licence GNU
GPL. Octave est un interpréteur de haut niveau, compatible la plupart du temps avec MATLAB et possédant la majeure partie
de ses fonctionnalités numériques. Dans ce chapitre, nous proposerons une introduction rapide à MATLAB et Octave. Le
but de ce chapitre est de fournir suffisamment d’informations pour pouvoir tester les méthodes numériques vues dans ce
polycopié. Il n’est ni un manuel de Octave/Matlab ni une initiation à la programmation.

1.1 Les environnements MATLAB et Octave

MATLAB et Octave sont des environnements intégrés pour le Calcul Scientifique et la visualisation. Ils sont écrits principale-
ment en langage C et C++. MATLAB est distribué par la société The MathWorks (voir le site www.mathworks.com). Son nom
vient de MATrix LABoratory, car il a été initialement développé pour le calcul matriciel. Octave, aussi connu sous le nom
de GNU Octave (voir le site www.octave.org), est un logiciel distribué gratuitement. Vous pouvez le redistribuer et/ou le
modifier selon les termes de la licence GNU General Public License (GPL) publiée par la Free Software Foundation.
Il existe des différences entre MATLAB et Octave, au niveau des environnements, des langages de programmation ou des
toolboxes (collections de fonctions dédiées à un usage spécifique). Cependant, leur niveau de compatibilité est suffisant pour
exécuter la plupart des programmes de ce cours indifféremment avec l’un ou l’autre. Quand ce n’est pas le cas – parce que les
commandes n’ont pas la même syntaxe, parce qu’elles fonctionnent différemment ou encore parce qu’elles n’existent pas
dans l’un des deux programmes – nous l’indiquons et expliquons comment procéder.
Nous utiliserons souvent dans la suite l’expression “commande MATLAB” : dans ce contexte, MATLAB doit être compris
comme le langage utilisé par les deux programmes MATLAB et Octave. De même que MATLAB a ses toolboxes, Octave
possède un vaste ensemble de fonctions disponibles à travers le projet Octave-forge. Ce dépôt de fonctions ne cesse de
s’enrichir dans tous les domaines. Certaines fonctions que nous utilisons dans ce polycopié ne font pas partie du noyau
d’Octave, toutefois, elles peuvent être téléchargées sur le site octave.sourceforge.net.

1.2 Installation(s) et version(s) en ligne

⋆ La documentation et les sources d’Octave peuvent être téléchargées à l’adresse https://www.gnu.org/software/
octave/.
La version en ligne d’Octave est disponible ici https://octave-online.net/.
⋆ L’université de Toulon propose aux étudiants la possibilité de le télécharger et de l’installer sur leur poste MATLAB.
Toutes les informations sont ici http://dsiun.univ-tln.fr/MATLAB.html
Par ailleurs, la version on line de MATLAB est disponible ici https://fr.mathworks.com/products/matlab-online.
html Les étudiants et enseignants de l’université de Toulon peuvent s’y connecter avec leurs paramètres universitaires.
Une fois qu’on a installé MATLAB ou Octave, on peut accéder à l’environnement de travail, caractérisé par le symbole d’invite
de commande >> . Il représente le prompt : cette marque visuelle indique que le logiciel est prêt à lire une commande. Il
suffit de saisir à la suite une instruction puis d’appuyer sur la touche «Entrée».

1.3 Premiers pas

Lorsqu’on démarre Octave, une nouvelle fenêtre va s’ouvrir, c’est la fenêtre principale qui contient trois onglets : l’onglet
“Fenêtre de commandes”, l’onglet “Éditeur” et l’onglet “Documentation”.

1.3.1 Fenêtre de commandes : mode interactif

L’onglet “Fenêtre de commandes” permet d’entrer directement des commandes et dès qu’on écrit une commande, Octave
l’exécute et renvoie instantanément le résultat. L’invite de commande se compose de deux chevrons (>>) et représente le
prompt : cette marque visuelle indique qu’Octave est prêt à lire une commande. Il suffit de saisir à la suite une instruction

3
Chapitre 1 Introduction à Octave/Matlab Mis à jour le Lundi 28 août 2023

puis d’appuyer sur la touche «Entrée». La console Octave fonctionne comme une simple calculatrice : on peut saisir une
expression dont la valeur est renvoyée dès qu’on presse la touche «Entrée». Voici un exemple de résolution d’un système
d’équations linéaires : 1
>> A = [2 1 0; -1 2 2; 0 1 4];
>> b = [1; 2; 3];
>> soln = A\b
soln =
0.25000
0.50000
0.62500

Ce mode interactif est très pratique pour rapidement tester des instructions et directement voir leurs résultats. Son utilisation
reste néanmoins limitée à des programmes de quelques instructions. En effet, devoir à chaque fois retaper toutes les
instructions s’avérera vite pénible.
Si on ferme Octave et qu’on le relance, comment faire en sorte que l’ordinateur se souvienne de ce que nous avons tapé ? On
ne peut pas sauvegarder directement ce qui se trouve dans la onglet “Fenêtre de commandes”, parce que cela comprendrait à
la fois les commandes tapées et les réponses du système. Il faut alors avoir préalablement écrit un fichier avec uniquement
les commandes qu’on a tapées et l’avoir enregistré sur l’ordinateur avec l’extension .m. Une fois cela fait, on demandera à
Octave de lire ce fichier et exécuter son contenu, instruction par instruction, comme si on les avait tapées l’une après l’autre
dans la Fenêtre de commandes. Ainsi plus tard on pourra ouvrir ce fichier et lancer Octave sans avoir à retaper toutes les
commandes. Passons alors à l’onglet “Éditeur”.

1.3.2 Éditeur : mode script

On voit qu’il n’y a rien dans cette nouvelle fenêtre (pas d’en-tête comme dans la “Fenêtre de commandes”). Ce qui veut
dire que ce fichier est uniquement pour les commandes : Octave n’interviendra pas avec ses réponses lorsque on écrira le
programme et ce tant que on ne le lui demandera pas. Ayant sauvé le programme dans un fichier avec l’extension .m, pour le
faire tourner et afficher les résultats dans la “Fenêtre de commandes” il suffira d’appuyer sur la touche «F5». Si on a fait une
faute de frappe, Octave le remarquera et demandera de corriger.
Maintenant qu’on a sauvé le programme, on est capable de le recharger.
Un fichier de script contient des instructions qui sont lues et exécutées séquentiellement par l’interpréteur d’Octave. Ce
sont obligatoirement des fichiers au format texte. Copier par exemple les lignes suivantes dans un fichier appelé first.m
A = [2 1 0; -1 2 2; 0 1 4];
b = [1; 2; 3];
soln = A\b

Appuyer sur la touche «F5», cliquer sur “Changer de répertoire” et regarder le résultat dans l’onglet “Fenêtre de commandes”. 2

1.4 Notions de base

1.4.1 Variables et affectation
Une variable peut être vue comme une boîte représentant un emplacement en mémoire qui permet de stocker une valeur et
à qui on a donné un nom afin de facilement l’identifier (boîte ← valeur) :

>> x=1 >> x=[2 5] >> x='c'

x = 1 x = x = c
2 5

L’affectation x=[2 5] crée une association entre le nom x et [2, 5]

le vecteur [2, 5] : la boîte de nom x contient le vecteur [2, 5].

³ 2 1 0 ´ µ x1 ¶ ³ 1 ´
1. Ces instructions calculent la solution du système linéaire −1 2 2 x2 = 2 . Noter l’usage des points-virgules à la fin de certaines instructions du
0 14 x3 3
fichier : ils permettent d’éviter que les résultats de ces instructions soit affiché à l’écran pendant l’exécution du script.
2. Sinon, si ce fichier se trouve dans le répertoire courant d’Octave, pour l’exécuter on peut juste taper son nom (sans l’extension) sur la ligne de
commande d’Octave : >> first
On peut aussi l’exécuter au moyen de la commande source qui prend en argument le nom du fichier ou son chemin d’accès (complet ou relatif au
répertoire courant). Par exemple : >> source("Bureau/TP1/first.m")

4 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 1.5 Commentaires

Il faut bien prendre garde au fait que l’instruction d’affectation (=) n’a pas la même signification que le symbole d’égalité
(=) en mathématiques (ceci explique pourquoi l’affectation de 1 à x, qu’en Octave s’écrit x = 1, en algorithmique se note
souvent x ← 1).
Une fois une variable initialisée, on peut modifier sa valeur en utilisant de nouveau l’opérateur d’affectation (=). La valeur
actuelle de la variable est remplacée par la nouvelle valeur qu’on lui affecte. Dans l’exemple précédent, on initialise une
variable à la valeur 1 et on remplace ensuite sa valeur par le vecteur [1, 2].
Il est très important de donner un nom clair et précis aux variables. Par exemple, avec des noms bien choisis, on comprend
tout de suite ce que calcule le code suivant :
base = 8
hauteur = 3
aire = base * hauteur / 2

Octave distingue les majuscules des minuscules. Ainsi mavariable, Mavariable et MAVARIABLE sont des variables diffé-
rentes.
Les noms de variables peuvent être non seulement des lettres, mais aussi des mots ; ils peuvent contenir des chiffres (à
condition toutefois de ne pas commencer par un chiffre), ainsi que certains caractères spéciaux comme le tiret bas «_» (appelé
underscore en anglais). Cependant, certains mots sont réservés :

ans Nom pour les résultats

eps Le plus petit nombre tel que 1+eps>1
inf ∞
NaN Not a number
i ou j i
pi π

>> 5/0
warning: division by zero
ans = Inf
>> 0/0
warning: division by zero
ans = NaN
>> 5*NaN % Most operations with NaN result in NaN
ans = NaN
>> NaN==NaN % Different NaN's are not equal!
ans = 0
>> eps
ans = 2.2204e-16

Si on écrit une instruction sans affectation, le résultat sera affecté à la variable ans.
>> [4,3]
ans =

4 3
>> 'Ciao'
ans = Ciao

Pour effacer la mémoire et désaffecter toutes les variables, utiliser la fonction clear all.

1.5 Commentaires
Le symbole % indique le début d’un commentaire : tous les caractères entre % et la fin de la ligne sont ignorés par
l’interpréteur.
⋆ Dans l’éditeur d’Octave, pour commenter plusieurs lignes en même temps, les sélectionner et appuyer sur les touches
«Ctrl+R». Pour dé-commenter plusieurs lignes en même temps, les sélectionner et appuyer sur les touches «Ctrl+Maj+R».
⋆ Dans l’éditeur de Matlab, pour commenter plusieurs lignes en même temps, les sélectionner et appuyer sur les touches
«Ctrl+R». Pour dé-commenter plusieurs lignes en même temps, les sélectionner et appuyer sur les touches «Ctrl+T».

1.6 Affichage
Lors de l’affectation d’une variable, le résultat de l’affectation sera affiché ; le symbole ; supprime cet affichage.

© 2022-2023 G. Faccanoni 5
Chapitre 1 Introduction à Octave/Matlab Mis à jour le Lundi 28 août 2023

>> a=[1,2]
a =

1 2

>> a=[4,3];
>> 'Ciao'
ans = Ciao

Pour afficher seulement le contenu d’une variable utiliser la fonction disp (en effet, si on écrit juste le nom de la variable, on
affichera aussi le nom de la variable)
>> a=[4,3];
>> disp(a)
4 3
>> a
a =

4 3
>> disp('Ciao')
Ciao

Pour nettoyer la fenêtre de commandes, utiliser la fonction clc.

1.7 Opérations arithmétiques

Dans Octave on a les opérations arithmétiques usuelles :

+ Addition
- Soustraction
* Multiplication
/ Division
ˆ Exponentiation

Quelques exemples :

>> a = 100 >> c = a-b >> a^b

a = 100 c = 83 ans = 1.0000e+34
>> b = 17 >> a/b
b = 17 ans = 5.8824

Les opérateurs arithmétiques possèdent chacun une priorité qui définit dans quel ordre les opérations sont effectuées. Par
exemple, lorsqu’on écrit 1 + 2 * 3, la multiplication va se faire avant l’addition. Le calcul qui sera effectué est donc 1 +
(2 * 3). Dans l’ordre, l’opérateur d’exponentiation est le premier exécuté, viennent ensuite les opérateurs *, /, // et %, et
enfin les opérateurs + et -.
Lorsqu’une expression contient plusieurs opérations de même priorité, ils sont évalués de gauche à droite. Ainsi, lorsqu’on
écrit 1 - 2 - 3, le calcul qui sera effectué est (1 - 2) - 3. En cas de doutes, vous pouvez toujours utiliser des parenthèses
pour rendre explicite l’ordre d’évaluation de vos expressions arithmétiques.
Il existe aussi les opérateurs augmentés (seulement dans Octave) :

a += b équivaut à a = a+b
a -= b équivaut à a = a-b
a *= b équivaut à a = a*b
a /= b équivaut à a = a/b
a ˆ = b équivaut à a = aˆ b

1.8 Division euclidienne

Lorsqu’on divise un nombre entier D (appelé dividende) par un autre nombre entier d (appelé diviseur), on obtient deux
résultats : un quotient q et un reste r , tels que D = qd + r (avec r < d ). La valeur q est le résultat de la division entière et la
valeur r celui du reste de cette division. Par exemple, si on divise 17 par 5, on obtient un quotient de 3 et un reste de 2 puisque
17 = 3 × 5 + 2. Ces deux opérateurs sont très utilisés dans plusieurs situations précises. Par exemple, pour déterminer si un
nombre entier est pair ou impair, il suffit de regarder le reste de la division entière par deux. Le nombre est pair s’il est nul et
est impair s’il vaut 1. Une autre situation où ces opérateurs sont utiles concerne les calculs de temps. Si on a un nombre de

6 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 1.9 Matrices

secondes et qu’on souhaite le décomposer en minutes et secondes, il suffit de faire la division par 60. Le quotient sera le
nombre de minutes et le reste le nombre de secondes restant. Par exemple, 175 secondes correspond à 175//60=2 minutes
et 175%60=55 secondes.
>> q=fix(9/4)
q = 2
>> % Reste de la division euclidienne de 9 par 4
>> r=rem(9,4)
r = 1
>>
>> r=mod(9,4) % 9 modulo 4
r = 1
>> q=fix(175/60)
q = 2
>> r=rem(175,60)
r = 55

1.9 Matrices
Pour définir une matrice on doit écrire ses éléments de la première à la dernière ligne, en utilisant le caractère ; pour séparer
les lignes (ou aller à la ligne). Notons que le symbole ; a deux fonctions : il supprime l’affichage d’un résultat intermédiaire
et il sépare les lignes d’une matrice. Par exemple, la commande
>> A = [ 1 2 3; 4 5 6]

ou la commande
>> A = [ 1 2 3
4 5 6]

donnent
A =
1 2 3
4 5 6

c’est-à-dire, une matrice 2 × 3 dont les éléments sont indiqués ci-dessus.

Un vecteur colonne est une matrice 1 × n, un vecteur ligne est une matrice n × 1 :
>> b = [1 2 3]
b =
1 2 3

>> b = [1; 2; 3]
b =
1
2
3

L’opérateur transposition s’obtient par la commande ' :

>> b = [1 2 3]'
b =
1
2
3

En Octave, les éléments d’une matrice sont indexés à partir de 1. Pour extraire les éléments d’une matrice on utilise la
commande A(i,j) où i et j sont la ligne et la colonne respectivement. On peut extraire un sous-vecteur en déclarant l’indice
i de début (inclus) et l’indice j de fin (inclus), séparés par deux-points v(i:j), ou encore un sous-vecteur en déclarant
l’indice i de début (inclus), le pas k et l’indice j de fin (inclus), séparés par des deux-points v(i:k:j). On peut même utiliser
un pas négatif. Cette opération est connue sous le nom de slicing (en anglais).
On peut combiner ces opérations pour extraire des sous-matrices :
A(2,3) % element A_{23}
A(:,3) % vecteur colonne [A_{13};...;A_{n3}]
A(1:4,3) % [A_{13};...A_{43}] premieres 4 lignes du vecteur colonne [A_{13};...A_{n3}]
A(1,:) % vecteur ligne [A_{11},...,A_{1n}]

© 2022-2023 G. Faccanoni 7
Chapitre 1 Introduction à Octave/Matlab Mis à jour le Lundi 28 août 2023

A(2,3:end) % [A_{23},...,A_{2n}] vecteur ligne

diag(A) % vecteur colonne [A_{11};...;A_{nn}] contenant la diagonale de A

Voici des exemples :

>> A = [8 1 6; 3 5 7; 4 9 2]
A =
8 1 6
3 5 7
4 9 2

>> A(2,3) % Element a la ligne 2 colonne 3

ans = 7
>> A(:,2) % Toutes les lignes, deuxieme colonne
ans =
1
5
9

>> A(2:3,2:3) % Sous-matrice 2 x 2

ans =
5 7
9 2

>> A(3:-1:1,:) % les lignes de la derniere a la premiere, toutes les colonnes

ans =
4 9 2
3 5 7
8 1 6

ATTENTION
Dans Octave les indices commencent à 1, ainsi A(1,:) indique la première ligne, A(2,:) la deuxième etc.

1.9.1 Concaténation de matrices

Nous pouvons générer une matrice par concaténation de deux ou plusieurs autres matrices :
⋆ Concaténation horizontale : A = [B C] ou horzcat(A,B)
⋆ Concaténation verticale : A = [B ; C] ou vertcat(A,B)
Dans le premier cas, on concatène côte à côte (horizontalement) les matrices B et C. Dans le second, on concatène vertica-
lement les matrices B et C. Attention aux dimensions qui doivent être cohérentes : dans le premier cas toutes les matrices
doivent avoir le même nombre de lignes, et dans le second cas le même nombre de colonnes.
Voici des exemples :
>> B = [1 2 3 ; 4 5 6]
B =

1 2 3
4 5 6

>> C = [7 8 9; 10 11 12]
C =

7 8 9
10 11 12

>> A = [B C]
A =

1 2 3 7 8 9
4 5 6 10 11 12

>> A = [B;C]
A =

8 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 1.9 Matrices

1 2 3
4 5 6
7 8 9
10 11 12

1.9.2 Matrices particulières

Construction de matrices particulières :
① La commande zeros(m,n) construit la matrice rectangulaire nulle O, i.e. celle dont tous les éléments ai j sont nuls
pour i = 1, . . . , m et j = 1, . . . , n.
La commande zeros(n) est un raccourci pour zeros(n,n).
② La commande ones(m,n) construit une matrice rectangulaire dont les éléments ai j sont égaux à 1 pour i = 1, . . . , m et
j = 1, . . . , n.
La commande ones(n) est un raccourci pour ones(n,n).
③ La commande eye(m,n) renvoie une matrice rectangulaire dont les éléments valent 0 exceptés ceux de la diagonale
principale qui valent 1.
La commande eye(n) (qui est un raccourci pour eye(n,n)) renvoie une matrice carrée de dimension n appelée
matrice identité et notée I.
④ La commande A=[] définit une matrice vide.
⑤ La commande diag(v) où v est un vecteur de n éléments renvoie une matrice carré de taille n dont les éléments
valent 0 exceptés ceux de la diagonale principale qui valent v.
⑥ Soit v un vecteur de n composantes. La commande diag(v) renvoie une matrice diagonale carrée de dimension n qui
contient v sur la diagonale principale ; la commande diag(v,1) renvoie une matrice carrée de dimension n + 1 qui
contient v sur la sur-diagonale principale etc.
Notons que diag(ones(1,4)) équivaut à eye(4).

>> Z=zeros(2,3) >> O=ones(3,2) >> E=eye(2,5) >> A=[] >> F=diag(v) >> G=diag(v,1)
Z = O = E = A = [](0x0) F = G =
0 0 0 1 1 Diagonal Matrix Diagonal Matrix 0 1 0
0 0 0 1 1 1 0 0 1 0 0 0
1 1 0 0 0 2 0 0 0 2
0 1 0 0 0 3 0
0 0 >> v=[1 2 3] 0 0 0
v = 3
1 2 3 0 0 0
0

Construction de vecteurs :
① x=[debut:pas:fin]
x fin −x debut x N −x 1
② x=linspace(debut,fin,N) (x a N points donc le pas h est N −1 = N −1 )

x = [-5 : 0.25 : 1] % x(k)= -5 + 0.25*(k-1), tant que k<=(fin-debut)/N

y = linspace(-5, 1, 25) % y(k)= -5 + h*(k-1), k=1,2,...,N avec h=(fin-debut)/(N-1)

Notons que la première instruction ne garantit pas que le dernier point soit pris, cela dépend du pas choisi (et des erreurs
d’arrondis) :
x = [0 : 0.4 : 1] % output : x=0.00000 0.40000 0.80000

Dans la première instructions on peut utiliser un pas négatif :

x = [1 : -0.4 : 0] % output : x=1.00000 0.60000 0.20000

Dimensions :
A = eye(3,4);
[r,c] = size(A) % r=nb de lignes et c=nb de colonnes de A
x = [0:10];
n = length(x) % n=nb d'elements de x
n = numel(x) % n=nb d'elements de x

© 2022-2023 G. Faccanoni 9
Chapitre 1 Introduction à Octave/Matlab Mis à jour le Lundi 28 août 2023

1.9.3 Opérations entre matrices

Opérations sur les matrices (lorsque les dimensions sont compatibles) :
⋆ Somme C = A + B, i.e. Ci j = Ai j + Bi j : C=A+B
⋆ Produit C = AB, i.e. Ci j = k=1 Ai k + Bk j : C=A*B NB il s’agit du produit matriciel !
P

⋆ Division à droite C = AB−1 : C=A/B

⋆ Division à gauche C = A−1 B : C=A\B (si B est un vecteur colonne alors C est un vecteur colonne solution du système
linéaire AC = B)
⋆ Élévation à la puissance C = AAA : C=A^3
⋆ Calcul du déterminant (si la matrice est carrée) : det(A)
⋆ Calcul de la matrice inverse (si la matrice est inversible) : inv(A)

>> A=[1 2 3; 4 5 6] >> B=ones(2,3) >> C=[1 2; 3 4; 5 >> D=eye(3,2) >> E=A(1:2,1:2)
A = B = 6] D = E =
1 2 3 1 1 1 C = Diagonal Matrix 1 2
4 5 6 1 1 1 1 2 1 0 4 5
3 4 0 1
5 6 0 0

>> A+B
ans =
2 3 4
5 6 7

>> A*C
ans =
22 28
49 64

>> A/B
ans =
1.00000 1.00000
2.50000 2.50000

>> b=[28;64]
b =
28
64

>> A\b
ans =
2.0000
4.0000
6.0000

>> A\B
ans =
-5.0000e-01 -5.0000e-01 -5.0000e-01
8.3267e-17 8.3267e-17 8.3267e-17
5.0000e-01 5.0000e-01 5.0000e-01

>> E^2
ans =
9 12
24 33

Quand on tente d’effectuer des opérations entre matrices de dimensions incompatibles on obtient un message d’erreur.

>> A+C
error: operator +: nonconformant arguments (op1 is 2x3, op2 is 3x2)

10 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 1.9 Matrices

1.9.4 ♥ Opérations pointées ♥

Quand il s’agit des opérations impliquant des multiplication (donc le produit mais aussi la division et l’élévation à la
puissance), la multiplication de deux matrices, avec les notations habituelles, ne signifie pas la multiplication élément par
élément mais la multiplication au sens mathématique du produit matriciel. C’est pour cela qu’Octave utilise deux opérateurs
distincts pour représenter la multiplication matricielle * et la multiplication élément par élément .* Le point placé avant
l’opérateur indique que l’opération est effectuée élément par élément. Les autres opérations de ce type sont la division à
droite et l’élévation à la puissance :
⋆ Produit Ci j = Ai j Bi j : C=A.*B NB il s’agit du produit d’Hadamard et non pas du produit matriciel
⋆ Division Ci j = Ai j /Bi j : C=A./B
⋆ Élévation à la puissance Ci j = A3i j : C=A.^3
Ces opérations sont à la base de la “vectorisation” car elles permettent le remplacement d’une boucle par une opération
matricielle pointée qui est généralement beaucoup plus performante.
>> A=[1 2 3; 4 5 6]
A =
1 2 3
4 5 6
>> B=[0 2 0; 0 0 1]
B =
0 2 0
0 0 1
>> A.*B
ans =
0 4 0
0 0 6
>> A*B
error: operator *: nonconformant arguments (op1 is 2x3, op2 is 2x3)

1.9.5 Opérateurs de comparaison et connecteurs logiques

Les opérateurs de comparaison renvoient 1 si la condition est vérifiée, 0 sinon. Ces opérateurs sont
On écrit < > <= >= == ~=
Ça signifie < > ≤ ≥ = ̸=
Bien distinguer l’instruction d’affectation = du symbole de comparaison ==.

ATTENTION
Les opérateurs de comparaison agissent élément par élément, ainsi lorsqu’on les applique à une matrice le résultat est une
matrice qui contient que des 0 ou 1 (parfois appelée “masque”). Par exemple
>> A = [1 2 3; 4 -5 6]; B = [7 8 9; 0 1 2];
>> A>B
ans =
0 0 0
1 0 1

On peut utiliser un masque pour remplacer seulement les éléments qui satisfont une conditions :
>> A = [1 -2 3; 4 -5 6];
>> A(A<0)=-100
A =
1 -100 3
4 -100 6

Les masques sont à la base de la “vectorisation” car permettent le remplacement d’un boucle avec des conditions par une
opération matricielle qui est généralement beaucoup plus performante.

E XEMPLE (M ASQUE )
Étant donné trois vecteurs :
⋆ hotels : une liste de noms d’hôtels
⋆ ratings : leurs notes dans une ville
⋆ cutoff : la note minimale
on souhaite afficher les noms des hôtels ayant une note supérieure ou égale au seuil.

© 2022-2023 G. Faccanoni 11
Chapitre 1 Introduction à Octave/Matlab Mis à jour le Lundi 28 août 2023

>> hotels =["CityLights";"SeaView";"MarketPlace";"ResortSpa";"Nightingale";"Clubadub";"SkylineView

";"MarinaBay";"ComfortFirst";"VillageValley"]; % vecteur colonne
>> ratings = [7.2;8.7;6.5;9.3;4.3;6.9;8.8;5.9;7.4;9.1]; % vecteur colonne
>> cutoff = 8;
>> good = hotels(ratings>=cutoff,:) % NB ":" pour selectionner toute la chaine de caracteres
good =

SeaView
ResortSpa
SkylineView
VillageValley

Pour combiner des conditions complexes (par exemple x > −2 et x 2 < 5), on peut combiner les opérateurs de comparaison
avec les connecteurs logiques :

On écrit Ça signifie
& et
| ou
~ non

Par exemple
>> (A > B) | (B > 5)
ans =
1 1 1
1 0 1

1.10 Fonctions
1.10.1 Fonctions prédéfinies
De très nombreuses fonctions sont déjà disponibles dans Octave/Matlab. Voici quelques exemples de fonction mathématique :
abs(-5)

sin(pi)
cos(pi)
tan(pi)

factorial(5) % output : 120

r=rem(11,3)

round(3.7) % output : 4
round(3.3) % output : 3
round(-3.7) % output : -4
round(-3.3) % output : -3

fix(3.7) % output : 3
fix(3.3) % output : 3
fix(-3.7) % output : -3
fix(-3.3) % output : -3

floor(3.7) % output : 3
floor(3.3) % output : 3
floor(-3.7) % output : -4
floor(-3.3) % output : -4

ceil(3.7) % output : 4
ceil(3.3) % output : 4
ceil(-3.7) % output : -3
ceil(-3.3) % output : -3

12 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 1.10 Fonctions

Remarque (Fonction vectorisée)

Une fonction vectorisée (“universelle” ou «ufunc», abréviation de universal function, est le terme exacte) est une fonction
qui peut s’appliquer terme à terme aux éléments d’un vecteur ou d’une matrice. Si f est une ufunc et si a = [a 1 , a 2 , . . . , a n ]
est un tableau, alors f (a) renvoie le tableau [ f (a 1 ), f (a 2 ), . . . , f (a n )]. En générale les fonctions prédéfinies sont vectorisées,
autrement dit si on applique la fonction à une matrice, elle renvoie une matrice de la même taille en ayant appliqué la
fonction à chaque élément.
x=[1:5]

sqrt(x) % output 1.0000 1.4142 1.7321 2.0000 2.2361

exp(x) % output 2.7183 7.3891 20.0855 54.5982 148.4132
log(x) % output 0.00000 0.69315 1.09861 1.38629 1.60944
log10(x) % output 0.00000 0.30103 0.47712 0.60206 0.69897
log2(x) % output 0.00000 1.00000 1.58496 2.00000 2.32193

ismember(2,x) % output 1
ismember(11,x) % output 0

Certaines fonctions sont spécifiques aux matrices :

A=[1 2; 3 4];
det(A)
inv(A)
trace(A)
size(A) % dimensions de A
x=[1 2 3];
length(x) % longueur d'un vecteur
numel(x) % nombre d'elements d'un vecteur

1.10.2 Définition d’une fonction

On peut définir nos propres fonctions au moyen de la commande function et les utiliser dans plusieurs scripts.
function [y1,...,yN] = myfunc(x1,...,xM)
instruction_1
instruction_2
...
[y1,...,yN] = ...
end

La structure type d’une fonction est la suivante :

⋆ on utilise la commande function dans laquelle on indique les arguments (x1,...,xM) et la valeur de retour
[y1,...,yN]
⋆ cette déclaration est suivie du corps de la définition qui est un bloc d’instructions à exécuter et se termine par le mot-clé
end ou endfunction
Quelques conseils :
⋆ Lorsque l’on définit une fonction dans un fichier, il est préférable de mettre un ; à la fin de chaque commande
constituant la fonction. Ainsi, on évitera l’affichage de résultats intermédiaires. Attention cependant à ne pas en mettre
sur la première ligne.
⋆ Lorsque l’on définit une fonction, il est préférable d’utiliser systématiquement les opérateurs terme à terme .* ./ et

.ˆ au lieu de * / et ˆ si l’on veut que cette fonction puisse s’appliquer à des scalaires, mais aussi à des tableaux.

ATTENTION
Pour éviter de surcharger une fonction déjà définie dans Matlab/Octave, prendre l’habitude d’appeler ses fonctions par
my...

Voir aussi https://fr.mathworks.com/help/matlab/matlab_prog/create-functions-in-files.html

Fichiers fonctions
Par convention, chaque définition de fonction est stockée dans un fichier séparé qui porte le nom de la fonction suivi de
l’extension .m Ces fichiers s’appellent des fichiers de fonction. Notez que c’est la même extension que les fichiers de scripts
mais, de plus, il faut absolument que le fichier s’appelle comme la fonction qu’il contient.
La structure type d’un fichier de fonction est la suivante :

© 2022-2023 G. Faccanoni 13
Chapitre 1 Introduction à Octave/Matlab Mis à jour le Lundi 28 août 2023

⋆ toute ligne commençant par un # ou un % est considérée comme un commentaire

⋆ les premières lignes du fichier sont des commentaires qui décrivent la syntaxe de la fonction. Ces lignes seront affichées
si on utilise la commande help myfunc
⋆ la fonction elle-même.

ATTENTION
Pour éviter toute confusion, utilisez le même nom pour le fichier de fonction et la première fonction du fichier. Matlab/Octave
associe votre programme au nom du fichier, pas au nom de la fonction. Les fichiers de script ne peuvent pas avoir le même
nom qu’une fonction du fichier.

Remarque (Sous-fonctions)
Un fichier de fonction peut en réalité contenir plusieurs fonctions déclarées au moyen de la commande function mais
seule la première définition est accessible depuis un script. Les autres définitions concernent des fonctions annexes (on dit
parfois des sous-fonctions) qui ne peuvent être utilisées que dans la définition de la fonction principale.

Voici un exemple qui prend en entrée un vecteur de valeurs et renvoie la moyenne et la déviation standard :

Fichier stat.m Script

% Cette fonction calcule la moyenne et la deviation values = [12.7, 45.4, 98.9, 26.6, 53.1];
% standard d'un vecteur x [average,stdeviation] = stat(values)
function [m,s] = stat(x)
n = numel(x);
m = sum(x)/n;
s = sqrt(sum((x-m).^2/n));
end

une fonction qui calcule l’aire d’un triangle en fonction des longueurs a, b et c des côtés grâce à la
À titre d’exemple, écrivonsp
formule de Héron : Aire = p(p − a)(p − b)(p − c) où p = (a + b + c)/2 est le demi-périmètre. On crée pour cela un fichier au
format texte appelé heron.m contenant les instructions suivantes
% Calcule l'aire s d'un triangle par la formule de Heron.
% a, b, c sont les longueurs des aretes.
function s = heron(a, b, c)
p = (a+b+c)/2;
s = sqrt(p*(p-a)*(p-b)*(p-c));
endfunction

La définition donnée ci-dessus peut être testée directement en chargeant le fichier heron.m avec la commande source et en
invoquant la fonction sur la ligne de commande. Par exemple :
>> source("Bureau/TP1/heron.m")
>> heron(3,5,4)
ans = 6

Fonctions locales dans un script - Matlab VS Octave

Les versions récentes de MATLAB permettent la création de fonctions directement dans un script. C’est très pratique pour de
petites fonctions “outils”, qui ne nécessitent pas d’être externalisées. En ajoutant des fonctions locales vous pouvez éviter de
créer et de gérer des fichiers de fonctions séparés. Cependant, la syntaxe est différente entre Matlab et Octave.
Matlab, depuis la version R2016b. Les fonctions locales peuvent apparaître dans n’importe quel ordre mais doivent être
placées à la fin du fichier, après le code du script. Voici un exemple :
x = 1:10;
n = numel(x);
avg = mymean(x,n)
med = mymedian(x,n)

function a = mymean(v,n)
% MYMEAN Local function that calculates mean of array.
a = sum(v)/n;
end

function m = mymedian(v,n)
% MYMEDIAN Local function that calculates median of array.

14 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 1.10 Fonctions

w = sort(v);
if rem(n,2) == 1
m = w((n + 1)/2);
else
m = (w(n/2) + w(n/2 + 1))/2;
end
end

Octave. Les fonctions locales peuvent apparaître dans n’importe quel ordre mais doivent être placées avant le code du
script et après l’instruction 1; Voici un exemple :

% Prevent Octave from thinking that this is a function file:

1;

function a = mymean(v,n)
% MYMEAN Local function that calculates mean of array.
a = sum(v)/n;
end

function m = mymedian(v,n)
% MYMEDIAN Local function that calculates median of array.
w = sort(v);
if rem(n,2) == 1
m = w((n + 1)/2);
else
m = (w(n/2) + w(n/2 + 1))/2;
end
end

x = 1:10;
n = numel(x);
avg = mymean(x,n)
med = mymedian(x,n)

1.10.3 Fonctions anonymes (lambda functions)

Les fonctions anonymes permettent de définir une fonction directement dans le script, à condition que la fonction se
compose d’une seule instruction. La syntaxe usuelle d’une fonction anonyme est
fun = @(arg1, arg2,...,argn) [expr] ; % crochets facultatifs si une seule expression

Nous utiliserons les fonctions anonymes surtout pour écrire directement la fonction dans le fichier de script sans créer un
fichier séparé en étant compatibles à la fois avec Matlab et avec Octave.
Une application courante des fonctions anonymes consiste à définir une expression mathématique.
f = @(x) 2*x ; % equivaut a definir f(x)=2x
f(2) % on evalue f(2) et on obtient 4

Cela permet entre autre de calculer rapidement une solution approchée d’une équation :
% on veut resoudre x=cos(x)
f = @(x) x.^2-2 ; % on pose f(x)=x^2-2
% fsolve( fct dont on cherche un zero , un point pas trop eloigne de la solution )
fsolve( f, 1 )
fsolve( f,-1 )

La contrainte d’utiliser une seule instruction n’empêche pas de calculer plusieurs résultats (car on peut renvoyer un vecteur)
ni d’écrire des boucles (grâce à l’utilisation des instructions pointées) ni des conditions (grâce à l’utilisation de masques, par
exemple pour la définition d’une fonction par morceaux, comme on verra à la page 22).

1.10.4 Arrayfun (listes de compréhension)

Les listes de compréhension sont une fonctionnalité puissante de Python qui permet de créer de nouvelles listes en appliquant
une certaine expression à chaque élément d’une séquence (comme une liste, un tuple ou une chaîne de caractères). La
syntaxe générale est
[ f(x) for x in xx ]

© 2022-2023 G. Faccanoni 15
Chapitre 1 Introduction à Octave/Matlab Mis à jour le Lundi 28 août 2023

où f (x) est l’évaluation d’une fonction f en chaque élément x de la séquence xx.

De manière analogue, en Matlab/Octave on peut utiliser arrayfun qui permet d’appliquer une fonction à chaque élément
d’un tableau. La syntaxe générale est
arrayfun( f, xx );

Exemples :

Python Octave/Matlab Bien-sur dans cet exemple simplifié on

peut vectoriser directement f :
f = lambda x : x**2 f = @(x) x^2;
f = @(x) x.^2;
xx = [1,3,10] xx = [1,3,10];
xx = [1,3,10];
[ f(x) for x in xx ] arrayfun( f, xx )
f(xx)

1.11 Graphes de fonctions R → R

Pour tracer le graphe d’une fonction f : [a, b] → R, il faut tout d’abord générer une liste de points x i où évaluer la fonction f ,
puis la liste des valeurs f (x i ) et enfin, avec la fonction plot, Octave reliera entre eux les points (x i , f (x i )) par des segments.
Plus les points sont nombreux, plus le graphe est proche du graphe de la fonction f .
Pour générer les points x i on peut utiliser
⋆ soit l’instruction linspace(a,b,n) qui construit la liste de n éléments

b−a
[a, a + h, a + 2h, . . . , b = a + nh] avec h =
n −1

x=linspace(1,5,5)
x =
1 2 3 4 5

⋆ soit l’instruction [a:h:b] qui construit la liste de n = E ( b−a

h ) + 1 éléments

[a, a + h, a + 2h, . . . , a + nh]

Dans ce cas, attention au dernier terme : b peut ne pas être pris en compte.
Voici un exemple avec une sinusoïde (en utilisant la fonction prédéfinie sin) :
x = linspace(-5,5,101); # x = [-5:0.1:5] with 101 elements
y = sin(x); # operation is broadcasted to all elements of the array
plot(x,y)

On obtient une courbe sur laquelle on peut zoomer, modifier les marges et sauvegarder dans différents formats.
Si la fonction n’est pas prédéfinie, il est bonne pratique de la définir pour qu’elle opère composante par composante lorsqu’on
lui passe un vecteur ou une matrice. Les opérations /, * et \^ agissant sur elle doivent être remplacées par les opérations
point correspondantes ./, .* et \.^ qui opèrent composante par composante.
Par exemple, on se propose de tracer la fonction

f : [−2; 2] → R
1
x 7→
1 + x2
Suivant la façon de définir la fonction, on pourra utiliser l’une des trois méthodes suivantes.
Méthode 1. En utilisant une fonction anonyme.
Dans un script ou dans la prompt on écrit les instructions suivantes :
f=@(x)[1./(1+x.^2)] % declaration de la fonction
x=[-2:0.5:2];
y=f(x); % evaluation en plusieurs points
plot(x,y) % affichage des points (x_i,y_i)

Méthode 2. En utilisant une fonction locale.

Dans un script on écrit les instructions suivantes (attention : syntaxe différente selon qu’on utilise Matlab ou
Octave) :

16 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 1.11 Graphes de fonctions R → R

Matlab Octave
x=[-2:0.5:2]; 1;
y=f(x); % evaluation en plusieurs points
plot(x,y) % affichage des points (x_i,y_i) function y=f(x)
y=1./(1+x.^2);
function y=f(x) end
y=1./(1+x.^2);
end x=[-2:0.5:2];
y=f(x); % evaluation en plusieurs points
plot(x,y) % affichage des points (x_i,y_i)

Méthode 3. En utilisant un fichier fonction.

On utilise deux fichiers :
3.1. Dans le fichier f.m on écrit la fonction informatique suivante
function y=f(x)
y=1./(1+x.^2);
end

3.2. Dans un script ou dans la prompt on écrit

x=[-2:0.5:2];
y=f(x); % evaluation en plusieurs points
plot(x,y) % affichage des points (x_i,y_i)

1.11.1 Plusieurs courbes sur le même repère

On peut tracer plusieurs courbes sur le même repère. Par défaut, MATLAB/Octave efface la figure avant chaque commande
de traçage plot. Vous pouvez tracer plusieurs lignes soit en écrivant plusieurs courbes dans une même instruction plot ou
à l’aide de la commande de mise en attente hold on. Tant que vous n’utilisez pas de suspension (hold off) ou que vous ne
fermez pas la fenêtre, tous les tracés apparaissent dans la fenêtre de la figure actuelle.
Par exemple, dans la figure suivante, on a tracé la même fonction : la courbe bleu correspond à la grille la plus grossière, la
courbe rouge correspond à la grille la plus fine :
a = linspace(-5,5,5); % a = [-5,-3,-1,1,3,5]
fa = sin(a);
b = linspace(-5,5,10); % b = [-5,-4,-3,...,5]
fb = sin(b);
c = linspace(-5,5,101); % c = [-5,-4.9,-4.8,...,5]
fc = sin(c);
plot(a,fa,b,fb,c,fc)
% la derniere ligne peut etre remplacee par
% hold on
% plot(a,fa)
% plot(b,fb)
% plot(c,fc)
% hold off

On peut spécifier la couleur et le type de trait, changer les étiquettes des axes, donner un titre, ajouter une grille, une légende
etc.
Par exemple, dans le code ci-dessous "r-" indique que la première courbe est à tracer en rouge (red) avec un trait continu, et
"g." que la deuxième est à tracer en vert (green) avec des points.
x = linspace(-5,5,101); # x = [-5,-4.9,-4.8,...,5] with
101 elements
y1 = sin(x); # operation is broadcasted to all elements
of the array
y2 = cos(x);
plot(x,y1,"r-",x,y2,"g.")
legend(['sinus';'cosinus'])
xlabel('abscisses')
ylabel('ordonnees')
title('Comparaison de sin(x) et cos(x)')
grid

Voir la table 1.1 et la documentation de Matlab pour connaître les autres options.

© 2022-2023 G. Faccanoni 17
Chapitre 1 Introduction à Octave/Matlab Mis à jour le Lundi 28 août 2023

linestyle= color= marker=

- solid line r red . points
- dashed line g green , pixel
: dotted line b blue o filled circles
-. dash-dot line c cyan v triangle down
m magenta ^ triangle up
y yellow > triangle right
w white < triangle left symbols
k black * star
+ plus
s square
p pentagon
x x
X x filled
d thin diamond
D diamond

TABLE 1.1 – Quelques options de plot

1.11.2 Plusieurs “fenêtres” graphiques

Avec figure() on génère une nouvelle fenêtres graphique :

x = [-pi:0.05*pi:pi];
figure(1)
plot(x, sin(x), 'r')
figure(2)
plot(x, cos(x), 'g')

1.11.3 Plusieurs repères dans la même fenêtre

La fonction subplot(x,y,z) subdivise la fenêtre sous forme d’une matrice de x lignes et y colonnes et chaque case est
numérotée, z étant le numéro de la case où afficher le graphe. La numérotation se fait de gauche à droite, puis de haut en bas,
en commençant par 1.

x = [-pi:0.05*pi:pi];

subplot(4,3,1)
plot(x, sin(x), 'r')

subplot(4,3,5)
plot(x, cos(x), 'g')

subplot(4,3,9)
plot(x, x.*x, 'b')

subplot(4,3,12)
plot(x, exp(-x.*x), 'm')

E XEMPLE
Dans le code suivant on voit comment tracer plusieurs courbes dans un même graphique (avec légende), plusieurs
graphe sur une même figure et plusieurs figures.

18 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 1.12 Structures itératives

figure(1)
x=[-2:0.5:2];

subplot(2,3,2)
plot(x,x,'r-',x,exp(x),'b*-')
legend(['y=x';'y=e^x'])
subplot(2,3,4)
plot(x,x.^2)
title('y=x^2')
subplot(2,3,5)
plot(x,x.^3)
xlabel('Axe x')
subplot(2,3,6)
plot(x,sqrt(x))

figure(2)
x=linspace(0.1,exp(2),100);
plot(x,log(x));

1.12 Structures itératives

Les structures de répétition se classent en deux catégories : les répétitions inconditionnelles pour lesquelles le bloc d’ins-
tructions est à répéter un nombre donné de fois et les répétitions conditionnelles pour lesquelles le bloc d’instructions est à
répéter autant de fois qu’une condition est vérifiée.

1.12.1 Répétition for

Lorsque l’on souhaite répéter un bloc d’instructions un nombre déterminé de fois, on peut utiliser un compteur actif,
c’est-à-dire une variable qui compte le nombre de répétitions et conditionne la sortie de la boucle.
La syntaxe de la commande for est schématiquement
for var = expression
instruction_1
instruction_2
end

expression peut être un vecteur ou une matrice. Par exemple, le code suivant calcule les 12 premières valeurs de la suite de
Fibonacci définie par la relation de récurrence u n = u n−1 + u n−2 avec pour valeurs initiales u 1 = u 2 = 1 :

n = 12; n = 12;
u = ones(1, n); # allocation u = [1,1];
for i = 3:n for i = 3:n
u(i) = u(i-1)+u(i-2); u = [ u , u(end)+u(end-1) ]; # concatenation
end end
disp(u) disp(u)

Dans les deux cas le résultat affiché est

1 1 2 3 5 8 13 21 34 55 89 144

Dans l’exemple suivant on calcul la somme des n premiers entiers (et on vérifie qu’on a bien n(n + 1)/2) :
n=100;
s=0;
for i=1:n
s += i;
end
s
n*(n+1)/2

Bien sur, dans ce cas il est préférable d’écrire

sum(1:100)

Il est possible d’imbriquer des boucles, c’est-à-dire que dans le bloc d’une boucle, on utilise une nouvelle boucle.

© 2022-2023 G. Faccanoni 19
Chapitre 1 Introduction à Octave/Matlab Mis à jour le Lundi 28 août 2023

for x = [10,20,30,40,50] % for x = 10:10:50

for y=[3,7]
disp(x+y)
end
end

Dans ce petit programme x vaut d’abord 10, y prend la valeur 3 puis la valeur 7 (le programme affiche donc d’abord 13, puis
17). Ensuite x = 20 et y vaut de nouveau 3 puis 7 (le programme affiche donc ensuite 23, puis 27).
Voir aussi https://fr.mathworks.com/help/matlab/ref/for.html

1.12.2 Boucle while : répétition conditionnelle

While est la traduction de “tant que. . .”. Concrètement, la boucle s’exécutera tant qu’une condition est remplie (donc tant
qu’elle renverra la valeur 1). Le constructeur while a la forme générale suivante :
i = ...
while condition(i)
instruction_1
instruction_2
i = ...
end

où condition représente des ensembles d’instructions dont la valeur est 1 ou 0. Tant que la condition condition a la
valeur 1, on exécute les instructions instruction_i.

ATTENTION
Si la condition ne devient jamais fausse, le bloc d’instructions est répété indéfiniment et le programme ne se termine pas.

Voici un exemple pour créer la liste 1, 12 , 13 , 14 :

£ ¤

nMax = 4;
n = 1;
a = [];
while n<=nMax
a=[a,1/n]; # Append element to list
n += 1; % si Matlab: n=n+1
end
disp(a)

Dans l’exemple suivant on calcul la somme des n premiers entiers tant que la somme ne dépasse pas 100 :
s=0;
n=0;
while s<100
n += 1; % si Matlab: n=n+1
s += n; % si Matlab: s=s+n
end
disp(n-1)
disp(s-n)
% En effet avec le dernier n on depasse 100

1.12.3 Vectorisation, i.e. optimisation des performances

La plupart du temps on manipule des vecteurs et des matrices. Les opérateurs et les fonctions élémentaires sont conçus pour
favoriser ce type de manipulation et, de manière plus générale, pour permettre la vectorisation des programmes. Certes, le
langage Octave contient des instructions conditionnelles, des boucles et la programmation récursive, mais la vectorisation
permet de limiter le recours à ces fonctionnalités qui ne sont jamais très efficaces dans le cas d’un langage interprété. Les
surcoûts d’interprétation peuvent être très pénalisants par rapport à ce que ferait un programme C ou FORTRAN compilé
lorsque l’on effectue des calculs numériques. Il faut donc veiller à réduire autant que possible le travail d’interprétation en
vectorisant les programmes.
Dans les exemples suivant, la fonction tic s’utilise avec la fonction toc pour mesurer le temps écoulé. La fonction tic
enregistre l’heure actuelle et la fonction toc utilise la valeur enregistrée pour calculer le temps écoulé.

E XEMPLE
Quasiment toutes les fonctions prédéfinies sont vectorisées. Voici un exemple avec la fonction sin.

20 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 1.13 Structure conditionnelle

Le code est significativement plus lent que

n = 100000;
xx = linspace(0,2*pi,n); n = 100000;
yy = zeros(numel(xx)); xx = linspace(0,2*pi,n);
tic tic
for i = 1:n yy = sin(xx);
yy(i) = sin(xx(i)); toc
end;
toc

E XEMPLE
1
Pour calculer 10000
P
n=1 n 2 , on peut utiliser les trois codes suivants, le deuxième étant significativement plus rapide :

n=1:10^7; n=1:10^7; n=1:10^7;

tic tic tic

s=0; s=sum(1./n.^2) s=(1./n)*(1./n)'
for i = n, toc toc
s+=1/i^2; % si Matlab: s=s+1/
i^2;
end;
s
toc

Un exemple de sorties obtenues :

⋆ avec le premier script : Elapsed time is 12.7138 seconds.
⋆ avec le deuxième script : Elapsed time is 0.116321 seconds.
⋆ avec le troisième script : Elapsed time is 0.098047 seconds.

1.13 Structure conditionnelle

Supposons vouloir calculer la valeur y = f (x) d’un nombre x selon la règle suivante :


 x si x ≤ −5,


 100 si − 5 < x ≤ 0,
f (x) =


 x2 si 0 < x < 10,

x −2 sinon.


On a besoin d’une instruction qui opère une disjonction de cas. En Octave il s’agit de l’instruction de choix introduite par le
mot-clé if. La syntaxe complète est la suivante :
if condition_1
instruction_1.1
instruction_1.2
...
elseif condition_2
instruction_2.1
instruction_2.2
...
...
else
instruction_n.1
instruction_n.2
...
end

où condition_1, condition_2. . . représentent des ensembles d’instructions dont la valeur est 1 ou 0 (on les obtient en gé-
néral en utilisant les opérateurs de comparaison). La première condition condition_i ayant la valeur 1 entraîne l’exécution
des instructions instruction_i.1, instruction_i.2. . . Si toutes les conditions sont 0, les instructions instruction_n
.1, instruction_n.2. . . sont exécutées. Les blocs elseif et else sont optionnels.
Pour l’exemple donné, la fonction (vectorisée) peut s’écrire comme suit :

© 2022-2023 G. Faccanoni 21
Chapitre 1 Introduction à Octave/Matlab Mis à jour le Lundi 28 août 2023

function y=f1(x)
n=numel(x);
for i=1:n
if x(i)<=-5
y(i)=x(i);
elseif x(i)<=0
y(i)=100;
elseif x(i)<10
y(i)=x(i)^2;
else
y(i)=x(i)-2;
end
end
end

La même fonction peut aussi s’écrire comme suit :

f2 = @(x) [ x.*(x<=-5) + 100.*(x>-5).*(x<=0) + (x.^2).*(x>0).*(x<10) + (x-2).*(x>=10) ];

Pour vérifier qu’on a bien la même fonction on compare le graphe des deux fonctions :
xx=[-7:0.1:12];
yy1=f1(xx);
yy2=f2(xx);
plot(xx,yy1,'r',xx,yy2,'b.')

Voici un exemple pour établir si un nombre est positif

a=-1.5;

if a < 0
disp('negative')
elseif a > 0
disp('positive')
else
disp(sign = 'zero')
end

Voici un exemple pour établir si un nombre est compris entre deux valeurs :
x = 10;
minVal = 2;
maxVal = 6;

if (x >= minVal) && (x <= maxVal) % equivaut a (x >= minVal) .* (x <= maxVal)
disp('Value within specified range.')
elseif (x > maxVal)
disp('Value exceeds maximum value.')
else
disp('Value is below minimum value.')
end

Voir aussi https://fr.mathworks.com/help/matlab/ref/if.html

E XEMPLE
Étant donné deux matrices d’entrée A et B, vérifier si on peut calculer le produit AB. Si c’est le cas, créez une matrice C
qui contient le produit AB, sinon, C doit contenir une chaîne de caractère contenant un message d’erreur.
function C = in_prod(A,B)
[rA,cA]=size(A);
[rB,cB]=size(B);
if cA==rB
C = A*B;
else
C = "Have you checked the inner dimensions?"
end
end

# TESTS
C=in_prod([1 2],[2;3])

22 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 1.13 Structure conditionnelle

C=in_prod(-5,100)
C=in_prod([1 2;3 4],[5;6])
C=in_prod([1 2 3; 4 5 6],[2 5;3 6])

© 2022-2023 G. Faccanoni 23
Chapitre 1 Introduction à Octave/Matlab Mis à jour le Lundi 28 août 2023

1.14 Exercices
Pensez à placer la commande clear all au début de vos scripts, de manière à nettoyer l’environnement de travail
(cela effacera toutes les variables en mémoire). Vous pouvez aussi utiliser la commande clc pour nettoyer la fenêtre de
commandes.

Pour chaque exercice, on écrira les instructions dans un fichier script nommé par exemple exo_1_3.m (sans espaces, ni
accents ni points sauf pour l’extension .m). On pourra bien-sûr utiliser le mode interactif pour simplement vérifier une
commande mais chaque exercice devra in fine être résolu dans un fichier script. Tous ces scripts devront se trouver dans
un dossier dont le nom ne contient ni espaces, ni accents, ni points.

Exercice 1.1
Soit les matrices
     
1 2 3 µ ¶ 1 1
2 −1 0
A = −1 0 1 B= u= x  v = 0 
−1 0 1
0 1 0 x2 1

1. Calculer tous les produits possibles à partir de A, B, u et v.

2. Calculer (A − I)7 et en extraire le coefficient en position (2, 3).
3. Calculer A−1 et la trace de A−1 (i.e. la somme des coefficients sur la diagonale).

Correction
Sans utiliser un module spécifique, il n’est pas possible de faire des calculs formels avec MATLAB/Octave, donc on ne peut
pas utiliser u sans donner une valeur numérique à x.
A=[1 2 3; -1 0 1; 0 1 0] % 3x3
B=[2 -1 0; -1 0 1] % 2x3
v=[1;0;1] % 3x1
% Les produits possibles sont
B*A % (2x3)*(3x3) -> 2x3
A*v % (3x3)*(3x1) -> 3x1
B*v % (2x3)*(3x1) -> 2x1
%
Id=eye(3)
D=(A-Id)^7 % c'est bien ^7 (produit matriciel) et non .^7
D(2,3) % -> 153
%
invA=A^(-1) % ou inv(A)
sum(diag(invA))

Exercice 1.2 (Vectorisation - if)

Considérons la fonction f : R → R définie par


 0 si x < 0,


x si 0 ≤ x < 1,
f (x) =
 2−x
 si 1 ≤ x ≤ 2,


0 si x > 2.


Écrire cette fonction de deux manières différentes et en afficher le graphe :

1. avec une instruction conditionnelle du type if ... elseif ... else (vectorisée) dans une fonction définie
avec function,
2. avec une instruction conditionnelle vectorisée dans une fonction anonyme.

Correction
Avec function La fonction vectorisée peut s’écrire comme suit :
function y=f1(x)
n=numel(x);
for i=1:n

24 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 1.14 Exercices

if x(i)<0
y(i)=0;
elseif x(i)<1
y(i)=x(i);
elseif x(i)<=2
y(i)=2-x(i);
else
y(i)=0;
end
end
end

Pour la tester, il fat d’abord la sauvegarder dans un fichier qui porte le même nom (ici le fichier devra s’appeler f1.m)
puis on écrira un autre fichier qui contient le script (qui se trouve dans le même dossier) pour la tester :
xx=[-1:0.1:3];
yy1=f1(xx);
plot(xx,yy1,'r')

Pour écrire un script qui contient directement cette fonction (sans utiliser deux fichiers séparés, l’un contenant la
fonction et l’autre le script), il faut décider s’il sera exécuté avec Matlab ou avec Octave :

Matlab Octave
xx=[-1:0.1:3]; 1;
yy1=f1(xx);
plot(xx,yy1,'r') function y=f1(x)
n=numel(x);
function y=f1(x) for i=1:n
n=numel(x); if x(i)<0
for i=1:n y(i)=0;
if x(i)<0 elseif x(i)<1
y(i)=0; y(i)=x(i);
elseif x(i)<1 elseif x(i)<=2
y(i)=x(i); y(i)=2-x(i);
elseif x(i)<=2 else
y(i)=2-x(i); y(i)=0;
else end
y(i)=0; end
end end
end
end xx=[-1:0.1:3];
yy1=f1(xx);
plot(xx,yy1,'r')

Anonyme La fonction anonyme vectorisée peut s’écrire comme suit :

f2 = @(x) [ 0*(x<0) + x.*(x>=0).*(x<1) + (2-x).*(x>=1).*(x<2) + 0.*(x>2) ];

Cette fonction peut être écrite directement avec le script et on a alors le même script en Octave ou en Matlab :
f2 = @(x) [ 0*(x<0) + x.*(x>=0).*(x<1) + (2-x).*(x>=1).*(x<2) + 0.*(x>2) ];

xx=[-1:0.1:3];
yy2=f2(xx);
plot(xx,yy2,'b.')

Exercice 1.3
Soit x un vecteur de n valeurs. On rappelle les définitions suivantes (pour l’estimation sur une population à partir d’un
échantillon) :
1X n
Moyenne arithmétique m = xi
n i =1
1 X n
Variance V = (x i − m)2
n − 1 i =1
p
Écart-type σ = V

© 2022-2023 G. Faccanoni 25
Chapitre 1 Introduction à Octave/Matlab Mis à jour le Lundi 28 août 2023

Médiane C’est la valeur qui divise l’échantillon en deux parties d’effectifs égaux. Soit y le vecteur qui contint les
composantes de x ordonné, alors

 y n/2 + y n/2+1

si n est pair,
médiane = 2
y si n est impair.
n/2+1

Écrire une fonction qui renvoie la moyenne, l’écart type et la médiane d’un vecteur donné en entré. Vérifier l’implémen-
tation sur un vecteur au choix en comparant avec les valeurs obtenues par les fonctions prédéfinies.

Correction
Dans cette correction on utilise deux fichiers séparés, l’un contenant la fonction et l’autre le script. On peut écrire la fonction
et le script dans un même fichier comme indiqué dans l’exercice précédent mais il faudra décider en amont si on exécutera le
script avec Matlab ou Octave.

’Fichier stat.m’ ’Script’

function [moy,et,med] = stat(x) clear all
n = length(x);
% MOYENNE = fonction predefinie mean(x) x = [0 0 8 1 1 2 2]
moy = sum(x)/n; [my_moy,my_std,my_mediane] = stat(x)
% ECART_TYPE = fonction predefinie std(x) octave_moy = mean(x)
et = sqrt(sum((x-moy).^2)/(n-1)); octave_std = std(x)
% MEDIANE = fonction predefinie median(x) octave_mediane = median(x)
y = sort(x);
if (rem(n,2) = = 0) % si n est un nombre pair x = [2 3 3 2]
med = (y(n/2)+y((n/2)+1))/2; [my_moy,my_std,my_mediane] = stat(x)
else octave_moy = mean(x)
med = y((n+1)/2); octave_std = std(x)
end octave_mediane = median(x)
end

Exercice 1.4
Un dispositif fournit un signal s(t ) = A sin(2πt + ϕ) avec A et ϕ inconnus. On mesure le signal à deux instants (en ms) :
s(0.5) = −1.76789123 et s(0.6) = −2.469394443. On posera α = A cos(ϕ) et β = A sin(ϕ).
1. Écrire et résoudre le système d’inconnues α et β. En déduire A et ϕ.
2. Tracer le signal et montrer qu’il passe par les points mesurés.

Correction
Rappel : sin(a + b) = sin(a) cos(b) + cos(a) sin(b) ainsi

s(t ) = A sin(2πt + ϕ) = A sin(2πt ) cos(ϕ) + cos(2πt ) sin(ϕ) = α sin(2πt ) + β cos(2πt ).

¡ ¢

On doit donc résoudre le système linéaire :

(
α sin(π) + β cos(π) = −1.76789123
α sin 56 π + β cos 65 π = −2.469394443
¡ ¢ ¡ ¢

qu’on écriture matricielle s’écrit

¡ 6 ¢ α = −1.76789123 ¡ 6 ¢ α = −1.76789123
µ ¶µ ¶ µ ¶ µ ¶µ ¶ µ ¶
sin(π) cos(π) ¡0 ¢ −1
i.e.
sin 56 π cos 5 π β sin 56 π cos 5 π β
¡ ¢
−2.469394443 −2.469394443

xx=[sin(2*pi*0.5) , cos(2*pi*0.5) ; sin(2*pi*0.6) , cos(2*pi*0.6)]\[-1.76789123;-2.469394443]

alpha=xx(1)
beta=xx(2)
π
p
Puisqu’on trouve α = β, alors cos(ϕ) = sin(ϕ), i.e. ϕ = 4 et A = 2α avec α ≈ 1.7679 :
s=@(t)[sqrt(2)*alpha*sin(2*pi*t+pi/4)];
tt=[0:0.01:1];
hold on

26 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 1.14 Exercices

plot(tt,s(tt),'r-') % la courbe
plot([0.5 0.6],[-1.76789123;-2.469394443],'b*') % les deux points
plot([0 0.5 0.5], [-1.76789123 -1.76789123 -3], 'g:') % trait pointille
plot([0 0.6 0.6], [-2.469394443 -2.469394443 -3], 'g:') % trait pointille
xlabel('t')
ylabel('s')
hold off
print("signal.jpg") % sauvagarde de la figure en .jpg

Exercice 1.5 (fsolve, plot)

Soit la fonction

f : [−10, 10] → R
x 3 cos(x) + x 2 − x + 1
x 7→ p
x 4 − 3x 2 + 127

1. Tracer le graphe de la fonction f en utilisant seulement les valeurs de f (x) lorsque la variable x prend successive-
ment les valeurs −10, −9.2, −8.4, . . . , 8.4, 9.2, 10 (i.e. avec un pas 0.8).
2. Apparemment, l’équation f (x) = 0 a une solution α voisine de 2. En utilisant le zoom, proposer une valeur
approchée de α.
3. Tracer de nouveau le graphe de f en faisant varier x avec un pas de 0.05. Ce nouveau graphe amène-t-il à corriger
la valeur de α proposée ?
4. Demander à Octave d’approcher α (fonction fsolve).

Correction
En utilisant un pas de 0.8 il semblerai que α = 1.89. En utilisant un pas de 0.05 il semblerai que α = 1.965. En utilisant la
fonction fsolve on trouve α = 1.9629.
clear all; clc;
f = @(x) [(x.^3 .*cos(x)+x.^2-x+1) ./ (x.^4-sqrt(3)*x.^2+127)] ;

subplot(1,2,1)
xx = [-10:0.8:10];
plot(xx,f(xx),'r-')
grid()

subplot(1,2,2)
xx = [-10:0.05:10];
plot(xx,f(xx),'r-')
grid()

fsolve(f,1.9)

© 2022-2023 G. Faccanoni 27
Chapitre 1 Introduction à Octave/Matlab Mis à jour le Lundi 28 août 2023

Exercice 1.6 (fsolve, plot)

On se propose ici d’utiliser Octave pour résoudre graphiquement des équations.

Considérons un cercle de rayon r . Si nous traçons un angle

A
t (mesuré en radians) à partir du centre du cercle, les deux B
rayons formant cet angle coupent le cercle en A et B . Nous
t

r
appelons a l’aire délimitée par la corde et l’arc AB (en bleu
r2
sur le dessin). Cette aire est donnée par a = (t − sin(t )).
2
Pour un cercle donné (c’est à dire un rayon donné), nous
choisissons une aire (partie en bleu) a. Quelle valeur de
l’angle t permet d’obtenir l’aire choisie ? Autrement dit,
connaissant a et r , nous voulons déterminer t solution de
l’équation
2a
= t − sin(t ).
r2

1. Résoudre graphiquement l’équation en traçant les courbes correspondant aux membres gauche et droit de
l’équation (pour a = 4 et r = 2). Quelle valeur de t est solution de l’équation ?
2. Comment faire pour obtenir une valeur plus précise du résultat ?

Correction Le graphe nous dit que la solution est entre 2 et 3. On peut

t=[0:pi/180:2*pi]; calculer une solution approchée comme suit :
rhs=t-sin(t);
a=4; fsolve( @(x) 2*a/(r^2)-(x-sin(x)) , 2.5)
r=2;
lhs=2*a/(r^2)*ones(size(t));
plot(t,rhs,t,lhs), grid

Jusqu’ici nous avons représenté des courbes engendrées par des équations cartésiennes, c’est-à-dire des fonctions
de la forme© y =ª f (x). Pour cela, nous générons d’abord un ensemble de valeurs { x i }i =0...N puis l’ensemble
de valeurs y i i =0...N avec y i = f (x i ). Il existe d’autres types de courbes comme par exemple les courbes
paramétrées (engendrées par des équations paramétriques). Les équations paramétriques de courbes planes sont
de la forme (
x = u(t ),
y = v(t ),

où u et v sont deux fonctions cartésiennes et le couple (x; y) représente les coordonnées d’un point de la courbe
paramétrée. La courbe engendrée par l’équation cartésienne y = f (x) est une courbe paramétrée car il suffit de
poser u(t ) = t et v(t ) = f (t ). Un type particulier de courbe paramétrique est constitué par les équations polaires
de courbes planes qui sont de la forme 3 (
x = r (t ) cos(t ),
y = r (t ) sin(t ).

Exercice 1.7 (Courbe parametrée)

Tracer la courbe papillon (t ∈ [0; 100]) :
( ¡ t ¢¢
x(t ) = sin(t ) e cos(t ) − 2 cos(4t ) − sin5 12
¡
¡ t ¢¢
y(t ) = cos(t ) e cos(t ) − 2 cos(4t ) − sin5 12
¡

3. r représente la distance de l’origine O du repère (O; i, j) au point (x; y) de la courbe, et t l’angle avec l’axe des abscisses.

28 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 1.14 Exercices

x = @(t) [ sin(t).*( exp(cos(t))-2*cos(4*t)-(sin(t/12)).^5 ) ];

Correction
y = @(t) [ cos(t).*( exp(cos(t))-2*cos(4*t)-(sin(t/12)).^5 ) ];
tt = [0:0.05:100];
plot(x(tt),y(tt))

Exercice 1.8 (Équation polaire)

Tracer la courbe papillon (t ∈ [0; 100]) :
(
x(t ) = r (t ) cos(t )
avec r (t ) = sin(7t ) − 1 − 3 cos(2t ).
y(t ) = r (t ) sin(t )

r = @(t) [sin(7*t)-1-3*cos(2*t)];
Correction
x = @(t) [ r(t).*cos(t) ];
y = @(t) [ r(t).*sin(t) ];
tt = [0:0.05:100];
plot(x(tt),y(tt))

Exercice 1.9 (Fractales)

La répétition de transformations permet de tracer des figures géométriques appelées fractales.
Considérons la suite définie par récurrence suivante :
µ ¶ µ ¶ µ ¶ µ ¶µ ¶ µ ¶
x 0 x m 11 m 12 x q1
= , = +
y 0 0 y n+1 m 21 m 22 y n q2 n

Pour chaque cas, tracer l’ensemble de points de coordonnées (x n , y n ).

µ ¶ µ ¶
1 1 −1 1
Exemple-1 on pose M = et q = ,
2 1 1 −3
Dragon de Heighway on choisit au hasard et de façon équiprobable α ∈ [0; 1[ puis on pose
µ ¶ µ ¶
1 1
M= α< M1 + α ≥ M2
2 2
et µ ¶ µ ¶
1 1
q= α< q1 + α ≥ q2
2 2
µ ¶ µ ¶ µ ¶ µ ¶
1 1 −1 1 −1 −1 0 1
avec M1 = , M2 = et q1 = , q2 = avec n = 0, . . . , 10000.
2 1 1 2 1 −1 0 0
Fougère de Barnsley on choisit au hasard et de façon équiprobable α ∈ [0; 1[ puis on pose
µ ¶ µ ¶ µ ¶ µ ¶
1 1 86 86 93 93
M= α< M1 + ≤α< M2 + ≤α< M3 + α ≥ M4
100 100 100 100 100 100
et µ ¶ µ ¶ µ ¶ µ ¶
1 1 86 86 93 93
q= α< q1 + ≤α< q2 + ≤α< q3 + α ≥ q4
100 100 100 100 100 100
µ ¶ µ ¶ µ ¶ µ ¶ µ ¶ µ ¶
1 0 0 1 85 4 1 20 −26 1 −15 28 0 0
avec M1 = , M2 = , M3 = , M4 = et q1 = , q2 = ,
100 0 16 100 −4 85 100 23 22 100 26 24 0 1.6
µ ¶
0
q3 = q2 , q4 = avec n = 0, . . . , 10000.
0.44
Arbre on choisit au hasard et de façon équiprobable α ∈ [0; 1[ puis on pose
µ ¶ µ ¶ µ ¶
1 1 2 2
M= α< M1 + ≤α< M2 + α ≥ M3
3 3 3 3

© 2022-2023 G. Faccanoni 29
Chapitre 1 Introduction à Octave/Matlab Mis à jour le Lundi 28 août 2023

et µ ¶ µ ¶ µ ¶
1 1 2 2
q= α< q1 + ≤α< q2 + α ≥ q3
3 3 3 3
µ ¶ µ ¶ µ ¶ µ ¶ µ 1 3 ¶
1 0 0 6 cos(ϑ) − sin(ϑ) 1 5 cos(ψ) −6 sin(ψ) 1 1 − 8 cos(ϑ)
avec M1 = , M2 = , M3 = et q1 = 2 2
, q = 51 3 ,
200 0 51 8 sin(ϑ) cos(ϑ) 8 5 sin(ψ) 6 cos(ψ) 0 2 200 − 8 sin(ϑ)
µ 1 5 ¶
2 − 16 cos(ψ)
q3 = 153 avec n = 0, . . . , 100, ϑ = − π8 et ψ = π5 .
5
1000 − 16 sin(ψ)
Source http://www.ac-grenoble.fr/maths/PM/Ressources/568/TP_fractales_Python.pdf

Correction

clc;
clear all;
cas=4;

% p et q vecteurs colonne, M matrice carree

transform = @(p,M,q) M*p+q;

% MAIN
n=1000*(cas==1)+10000*(cas==2)+10000*(cas==3)+100*(cas==4);
A = zeros(n,2);
for i=2:n
if cas==1 % Exemple1
M=[1,-1;1,1]/2;
q=[1;-3];
elseif cas==2 % Dragon
M1=[1,-1;1,1]/2;
q1=[0;0];
M2=[-1,-1;1,-1]/2;
q2=[1;0];
alpha=rand();
M=(alpha<0.5)*M1+(alpha>=0.5)*M2;
q=(alpha<0.5)*q1+(alpha>=0.5)*q2;
elseif cas==3 % Fougere
M1=[0,0;0,0.16];
q1=[0;0];
M2=[0.85,0.04;-0.04,0.85];
q2=[0;1.6];
M3=[0.2,-0.26;0.23,0.22];
q3=[0;1.6];
M4=[-0.15,0.28;0.26,0.24];
q4=[0;0.44];
alpha=rand();
M=(alpha<0.01)*M1+(alpha>=0.01)*(alpha<0.86)*M2+(alpha>=0.86)*(alpha<0.93)*M3+(alpha>=0.93)*M4;
q=(alpha<0.01)*q1+(alpha>=0.01)*(alpha<0.86)*q2+(alpha>=0.86)*(alpha<0.93)*q3+(alpha>=0.93)*q4;
else % arbre
M1=[0,0;0,51/200];
q1=[0.5;0];
theta=-pi/8;
M2=6/8*[cos(theta),-sin(theta);sin(theta),cos(theta)];
q2=[1/2-3/8*cos(theta);51/200-3/8*sin(theta)];
psi=pi/5;
M3=1/8*[5*cos(psi),-6*sin(psi);5*sin(psi),6*cos(psi)];
q3=[1/2-5/16*cos(psi);153/1000-5/16*sin(psi)];
alpha=rand();
M=(alpha<1/3)*M1+(alpha>=1/3)*(alpha<2/3)*M2+(alpha>=2/3)*(alpha<0.93)*M3;
q=(alpha<1/3)*q1+(alpha>=1/3)*(alpha<2/3)*q2+(alpha>=2/3)*(alpha<0.93)*q3;
end
A(i,:)=transform(A(i-1,:)',M,q);
end
plot(A(:,1),A(:,2),'o','MarkerSize',8)

30 © 2022-2023 G. Faccanoni
Chapitre 2
Traitement mathématique des images numériques
Dans cette partie nous allons nous intéresser à la manipulation d’images numériques.
Pour commencer, nous allons considérer une matrice de 0 et 1 comme une image constituée de “carreaux” blancs si la valeur
est 0 et noirs si la valeur est 1. Nous allons alors appliquer quelque transformation élémentaire sur la matrice initiale et
visualiser sur l’image associée l’effet de chaque transformation.

Attention
Pensez à placer la commande clear all au début de vos scripts, de manière à nettoyer l’environnement de travail
(cela effacera toutes les variables en mémoire). Vous pouvez aussi utiliser la commande clc pour nettoyer la fenêtre de
commandes.

Pour chaque exemple ou exercice, on écrira les instructions dans un fichier script nommé par exemple exo_1_3.m
(sans espaces, ni accents ni points sauf pour l’extension .m). On pourra bien-sûr utiliser le mode interactif pour
simplement vérifier une commande mais chaque exercice devra in fine être résolu dans un fichier script. Tous ces
scripts devront se trouver dans un dossier dont le nom ne contient ni espaces, ni accents, ni points.

Exercice 2.1 (Multiplication matricielle appliquée)

On modélise une image en noir et blanc formée de 25 pixels par une matrice de 5 lignes et 5 colonnes, dans laquelle 0
correspond à un pixel blanc et 1 à un pixel noir. L’image à modéliser est la suivante :

1. Donner la matrice M associée à l’image.

2. Soient
0 0 0 0 1 0 1 0 0 0
   
0 0 0 1 0 0 0 1 0 0
A= B=
   
0 0 1 0 0 0 0 0 1 0
0 1 0 0 0 0 0 0 0 1
1 0 0 0 0 1 0 0 0 0
Calculer le produit AM. Quel est l’effet de la matrice A sur l’image ? Quel est l’effet si on fait le produit MA sur
l’image ? Quelles instructions permettent d’obtenir ces résultats sans utiliser la multiplication matricielle ?
Calculer le produit BM. Quel est l’effet de la matrice B sur l’image ? Quel est l’effet si on fait le produit MB sur
l’image ? Quelles instructions permettent d’obtenir ces résultats sans utiliser la multiplication matricielle ?
3. Quelle image obtient-on en faisant le produit AMB ? Quelles instructions permettent d’obtenir ces résultats sans
utiliser la multiplication matricielle ?
4. Quelle image obtient-on en faisant le produit AMA ? Quelles instructions permettent d’obtenir ces résultats
sans utiliser la multiplication matricielle ?
5. Quel produit matriciel peut-on faire pour obtenir la figure suivante ? Quelles instructions permettent d’obtenir
ces résultats sans utiliser la multiplication matricielle ?

31
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

6. Calculer le produit AMT . Quelles instructions permettent d’obtenir ces résultats sans utiliser la multiplication
matricielle ?

Correction
1.
0 1 1 1 0
 
0 1 0 0 0
M=
 
0 1 1 0 0

0 1 0 0 0
0 1 0 0 0

2. Le produit AM correspond à inverser l’ordre des lignes de la matrice M : l’image est alors symétrique par rapport à la
troisième ligne.
Le produit MA correspond à inverser l’ordre des colonnes de la matrice M : l’image est alors symétrique par rapport à
la troisième colonne.
Le produit BM correspond à translater les lignes de la matrice M d’un rang vers le haut (la première ligne passant en
cinquième ligne) : l’image est alors translatée d’un rang vers le haut (la première ligne passant en cinquième ligne).
Le produit MB correspond à translater les colonnes de la matrice M d’un rang vers la droite (la cinquième colonne
passant en première colonne) : l’image est alors translatée d’un rang vers la droite (la cinquième colonne passant en
première colonne).

0 1 0 0 0 0 1 1 1 0 0 1 0 0 0 0 0 1 1 1
       
0 1 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0
AM =  MA =  BM =  MB = 
       
0 1 1 0 0
 0 0 1 1 0 0 1 0 0 0 0 0 1 1 0

0 1 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0
0 1 1 1 0 0 0 0 1 0 0 1 1 1 0 0 0 1 0 0

M AM MA BM MB

3. Le produit AMB = (AM)B correspond par exemple à une symétrie horizontale par rapport à la troisième ligne suivie
d’une translation d’un rang vers la droite. On peut aussi l’écrire comme AMB = A(MB) qui correspond à une translation
d’un rang vers la droite suivie d’une symétrie horizontale par rapport à la troisième ligne. Dans tous les cas on obtient
l’image suivante :

M → AM → (AM)B

M → MB → A(MB)

4. On peut par exemple calculer AMAB.

M → AM → (AM)A → (AMA)B

Ou encore calculer AMBABB.

AMB → AMBA → AMBABB

32 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023

5. Le produit AMT correspond à une rotation antihoraire, elle équivaut à l’instruction rot90(M).

AMT

Par multiplications :

M=[0 1 1 1 0 A=eye(5); A*M

0 1 0 0 0 A=A(:,5:-1:1) M*A
0 1 1 0 0 % B*M
0 1 0 0 0 B=diag(ones(4,1),1); M*B
0 1 0 0 0] B(5,1)=1; A*M*B
B A*M*A*B
A*M*B*A*B*B

Par transformations :

M=[0 1 1 1 0 AM = M(end:-1:1,:) % octave flipud(M), matlab flip(

0 1 0 0 0 M)
0 1 1 0 0 MA = M(:,end:-1:1) % octave fliplr(M), matlab flip(
0 1 0 0 0 M,2)
0 1 0 0 0] BM = [ M(2:end,:) ; M(1,:) ]
% MB = [ M(:,end), M(:,1:end-1) ]
A=eye(5); AMB = [ AM(:,end), AM(:,1:end-1) ]
A=A(:,5:-1:1)
% AMA = AM(:,end:-1:1);
B=diag(ones(4,1),1); AMAB = [ AMA(:,end), AMA(:,1:end-1) ]
B(5,1)=1;
B AMt = M'(end:-1:1,:) % rot90(M)

© 2022-2023 G. Faccanoni 33
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

2.1 Introduction : lecture, affichage, sauvegarde

Un peu de nomenclature :
⋆ Acquisition d’une image (= numérisation) : dans la figure ci dessous on a à gauche une image réelle et à droite sa version
numérisée. La numérisation d’une image réelle comporte une perte d’information due d’une part à l’échantillonnage
(procédé de discrétisation spatiale d’une image consistant à projeter sur une grille régulière, i.e. en un nombre fini de
points, une image analogique continue en lui associant une valeur unique), d’autre part à la quantification (nombre
finis de nuances = limitation du nombre de valeurs différentes que peut prendre un point).

⋆ Les pixels d’une image : une image numérique en niveaux de gris (grayscale image en anglais) est un tableau de valeurs
A. Chaque case de ce tableau, qui stocke une valeur a i j , se nomme un pixel (PICTure ELement).
En notant n le nombre de lignes et p le nombre de colonnes du tableau, on manipule ainsi un tableau de n × p pixels.
Les valeurs des pixels sont enregistrées dans l’ordinateur ou l’appareil photo numérique sous forme de nombres entiers
entre 0 et 255, ce qui fait 256 valeurs possibles pour chaque pixel. La valeur 0 correspond au noir et la valeur 255
correspond au blanc. Les valeurs intermédiaires correspondent à des niveaux de gris allant du noir au blanc.
On peut définir une fonction comme suit :

g : ‚1 : nƒ × ‚1 : pƒ → ‚0 : 255ƒ
(i , j ) 7→ g (i , j ) = a i j

⋆ La taille d’une image : la taille d’une image est le nombre de pixel. Les dimensions d’une image sont la largeur (=nombre
de colonnes du tableau) et la hauteur (=nombre de lignes du tableau).
⋆ La résolution d’une image : la résolution d’une image est le nombre de pixel par unité de longueur. En générale, on
utilise des “pixel par pouce” ou “point par pouce” (ppp), on anglais on dit dot per inch (dpi) : n dpi = n pixel pour un
puce = n pixel pour 2.54 cm.

E XEMPLE
On veut déterminer les dimensions d’une image obtenue à partir d’un scanner pour une page A4 à la résolution de 300
dpi.
⋆ Le format A4 est un rectangle de 21 cm × 29.7 cm = 8.25 inch × 11.7 inch.
⋆ La largeur de l’image est donc n = 300 dpi × 8.25 inch = 2 475 pixel.
⋆ La hauteur de l’image est donc p = 300 dpi × 11.7 inch = 3 510 pixel.
⋆ La taille de l’image est 2 475 × 3 510 = 8 700 000 pixel.

Matrice → Image
Considérons une matrice carrée avec n = p = 8, ce qui représente 23 × 23 = 26 = 64 éléments.

N = 8; A =
A = eye(N); -20 0 20 40 60 80 100 120
for n = 1:N 0 20 40 60 80 100 120 140
A(:,n) = 20*(n-2:n+N-3); 20 40 60 80 100 120 140 160
end 40 60 80 100 120 140 160 180
60 80 100 120 140 160 180 200
80 100 120 140 160 180 200 220
100 120 140 160 180 200 220 240
120 140 160 180 200 220 240 260

Pour voir le tableau A comme une image de 64 pixels en niveaux de gris,

34 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.1 Introduction : lecture, affichage, sauvegarde

⋆ d’abord on transforme la matrice avec uint8 (toutes les valeurs inférieurs à 0 sont mises à 0, toutes les valeurs
supérieures à 255 sont mises à 255 et les valeurs non entières sont arrondies)
⋆ on indique la colormap (ici échelle de gris de 0 à 255) pour que l’image s’affiche avec des niveaux de gris (si M i j = 255
est affiché blanc et M i j = 0 noir, les valeurs intermédiaires en gris), puis on affiche l’image :

M = uint8(A) M =
% si A(i,j)<0 alors M(i,j)=0, 0 0 20 40 60 80 100 120
% si A(i,j)>255 alors M(i,j)=255 0 20 40 60 80 100 120 140
20 40 60 80 100 120 140 160
% methode 1 40 60 80 100 120 140 160 180
%colormap(gray(256)); 60 80 100 120 140 160 180 200
%image(M) 80 100 120 140 160 180 200 220
100 120 140 160 180 200 220 240
% methode 2 120 140 160 180 200 220 240 255
imshow(M,gray(256))

% colorbar % pour voir la table de couleurs

utilisees (et les valeurs presentes dans l
image)

2 0.8

0.6
4

0.4

0.2

0
2 4 6 8

Image → Matrice
Dans ce cas, nous ne construisons pas la matrice à la main, mais nous allons lire un fichier image et l’importer comme une
matrice dans Octave. Pour transformer une image en une matrice il suffit d’indiquer dans un script : 1
A = imread('lena512.bmp');

ATTENTION
L’image doit se trouver dans le même dossier que les fichier script. Ce dossier et le fichier script ne doivent pas contenir ni
d’espaces, ni d’accents, ni de points ou autres caractères spéciaux.
Pour connaître la dimension de la matrice (i.e. de l’image) il suffira d’écrire
[row,col] = size(A) % = [hateur,largeur] de l'image

On a une matrice carrée de taille n = p = 512, ce qui représente 512 × 512 = 218 = 262 144 pixels. 2
Pour connaître l’intervalle des valeurs de la matrice, i.e. min A i j et max A i j , on écrira
1≤i ≤row 1≤i ≤row
1≤ j ≤col 1≤ j ≤col

min(A(:)), max(A(:))

Les niveaux de gris de notre image sont compris entre 25 et 245.

Par défaut, avec la fonction imread, Octave transforme un fichier image en niveaux de gris (ici de type .bmp) en une matrice
dont chaque coefficient est de type uint8 (codage sur 8 bits non signés). Cela signifie que ses éléments (qui représentent les
niveaux de gris) sont dans l’intervalle d’entiers ‚0 − 255ƒ). Pour pouvoir utiliser toutes les fonctions Octave sur ces données, il
est préférable de les convertir en réels (commande double) puis, après traitement, les reconvertir en entier codé sur 8 bits
non signé avant de sauvegarder en format raster (commande uint8). 3

1. Cette image, “Lena”, est une image digitale fétiche des chercheurs en traitement d’images. Elle a été scannée en 1973 dans un exemplaire de Playboy
et elle est toujours utilisée pour vérifier la validité des algorithmes de traitement ou de compression d’images. On la trouve sur le site de l’University of
Southern California http://sipi.usc.edu/database/.
2. Les appareils photos numériques peuvent enregistrer des images beaucoup plus grandes, avec plusieurs millions de pixels.
3. Octave peut lire des images codées sur 8, 16, 24 ou 32 bits. Mais le stockage et l’affichage de ces données ne peut être fait qu’avec trois types de
variables :
⋆ le type uint8 (entier non signé de 8 bits) de plage [0; 255] ;
⋆ le type uint16 (entier non signé de 16 bits) de plage [0; 65535] ;
⋆ le type double (réel 64 bits) de plage [0; 1].

© 2022-2023 G. Faccanoni 35
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

Matrice/Image → affichage et/ou sauvegarde

On peut voir la matrice comme une image avec :

imshow(A); % imshow(uint8(A)) si on veut s'assurer

d'avoir des uint8

On peut sauvegarder l’image avec :

imwrite(uint8(A),'monimage.jpg','jpg'); F IGURE 2.1 – Léna (original)

http://www.lenna.org

Exercice 2.2
Afficher les images associées aux matrices suivantes :

0 0 0 255 255
 
 
 0 0 255 255 255 0 50 100 150 200 250
A= B = 0
 
 0 255 255 255 0 
 50 100 150 200 250
255 255 255 0 0  0 50 100 150 200 250
255 255 0 0 0

Correction
Les valeurs 255 sont affichées en blanc, les valeurs 0 en noir, les valeurs intermédiaires en gris.
% Dans OCTAVE on peut ecrire une seule instruction
A=255*(eye(5)+diag(ones(1,4),1)+diag(ones(1,4),-1))(:,end
:-1:1)
% Dans MATLAB il faut la decouper comme suit
% A=255*(eye(5)+diag(ones(1,4),1)+diag(ones(1,4),-1))
% A=A(:,end:-1:1)

B=[0:50:250].*ones(3,1)

subplot(1,2,1)
imshow(A,gray(256));
title("A")
subplot(1,2,2)
% noter la difference avec imshow(B)
imshow(B,gray(256));
title("B")

Dans les prochaines sections nous allons considérer deux familles de transformations :
⋆ les transformations qui déplacent/éliminent/ajoutent des lignes/colonnes de la matrice ;
⋆ les transformations qui modifient la valeur des pixels.
Dans le deuxième cas, il convient de transformer d’abord les valeurs de type uint en double car certaines opérations ne
sont pas définies pour les uint. On retransformera en uint avant de sauvegarder l’image.

Remarque (uint8)
Lorsqu’on effectue des calculs avec des uint8, les valeurs sont toujours des entiers dans [0; 255] (toutes les valeurs inférieures
à 0 sont mises à 0, toutes les valeurs supérieures à 255 sont mises à 255 et les valeurs non entières sont arrondies). Dans le
code ci-dessous, on a a = 155 mais a est de type uint8 donc, lorsqu’on calcule a 2 , Octave plafonne à 255 ; lorsqu’on calcule
a/2, Octave l’arrondi à 78 et lorsqu’on calcule −a, Octave plafonne à 0 :
A = imread('lena512.bmp');
a = A(10,10);

% On verifie que par default imread donne une matrice de uint8

disp(a) % output 155
disp(class(a)) % output uint8
disp(a^2) % output 255 au lieu de 155^2=24025
disp(a/2) % output 78 au lieu de 155/2=77.5

36 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.1 Introduction : lecture, affichage, sauvegarde

disp(-a) % output 0 au lieu de -155

% On transforme en float et on verifie que les resultats ne sont plus tronques ni arrondis
a = double(a);
disp(a) % output 155
disp(class(a)) % output double
disp(a^2) % output 155^2=24025
disp(a/2) % output 155/2=77.5
disp(-a) % output -155

Les fonctions Octave les plus utiles pour gérer les images sont les suivantes :
⋆ imread : lit une image à partir d’un fichier (formats standards) ;
⋆ image ou imagesc ou imshow : affiche une image (objet graphique Image) ;
⋆ imwrite : écrit une image dans un fichier (formats standards) ;
⋆ imfinfo : extrait des informations d’un fichier (formats standards) ;
⋆ print : exporte une image (formats standards).

Remarque (image, imagesc ou imshow ?)

Soit
 
0 2 4 6
A= 8 10 12 14
16 18 20 22
Elle contient 3 lignes et 4 colonnes et les valeurs vont de 0 à 22. Analysons le code suivant :

clear all
image(A)
close all 0.5 1
1 0.8
1.5 0.6
2 0.4
2.5
A=[0 2 4 6; 8 10 12 14; 16 18 20 22]; 3
3.5
0.2
0
1 2 3 4
minA=min(A(:)) image(A,’CDataMapping’,’scaled’)
maxA=max(A(:)) 0.5
1 20
1.5 15
2 10
2.5 5
3
n=5 3.5
1 2 3 4
0

imagesc(A)
0.5 20
subplot(n,1,1) 1
1.5 15
2 10
2.5
image(A), colormap(gray(256)), colorbar 3 5
3.5 0
title("image(A)"); 1 2
imshow(A,[])
3 4

20
15
subplot(n,1,2) 10
5
image(A,'CDataMapping','scaled'), colormap(gray 0

(256)), colorbar imshow(A)

1
title("image(A,'CDataMapping','scaled')"); 0.8
0.6
0.4
0.2
0
subplot(n,1,3)
imagesc(A), colormap(gray(256)), colorbar
title("imagesc(A)");

subplot(n,1,4)
imshow(A,[]), colorbar
title("imshow(A,[])");

subplot(n,1,5)
imshow(A), colorbar
title("imshow(A)");

⋆ La fonction image affiche une image noire : en effet les valeurs de A vont de 0 à 22 tandis que l’échelle de gris va de 0 à
255. Pour étirer l’échelle il faut le demander explicitement avec les mots clé ’CDataMapping’,’scaled’.
⋆ La fonction imagesc affiche une image en “étirant” l’échelle : elle étale les valeurs pour que la plus petite valeur de A
correspond au 0 de l’échelle de gris et la plus grande valeur de A au 255 dans l’échelle de gris.
⋆ La fonction image choisi automatiquement la colormap et affiche les cellules comme des carrés mais a le même
problème que image. Pour qu’elle étale les valeurs pour que la plus petite valeur de A correspond au 0 de l’échelle de
gris et la plus grande valeur de A au 255 dans l’échelle de gris il faut ajouter [].
Voici un autre exemple avec une image bitmap et des options pour avoir le bon comportement :

© 2022-2023 G. Faccanoni 37
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

clear all, close all

image(A)
A=imread('lena512.bmp');
minA=min(A(:)) 1
0.8
maxA=max(A(:)) 0.6
0.4
0.2
n=5 0

image(A,’CDataMapping’,’scaled’)
subplot(n,1,1)
image(A), colormap(gray(256)), axis image, axis off 200
, colorbar 150
title("image(A)"); 100
50

subplot(n,1,2)
imagesc(A)
image(A,'CDataMapping','scaled'), colormap(gray
(256)), axis image, axis off, colorbar 200
title("image(A,'CDataMapping','scaled')"); 150
100
50
subplot(n,1,3)
imagesc(A), colormap(gray(256)), axis image, axis imshow(A,[])
off, colorbar
title("imagesc(A)"); 200
150
100
subplot(n,1,4) 50
imshow(A,[]), colorbar
title("imshow(A,[])"); imshow(A)
250
subplot(n,1,5) 200
150
imshow(A), colorbar 100
title("imshow(A)"); 50
0

38 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.2 Manipulations élémentaires : déplacements des pixels

2.2 Manipulations élémentaires : déplacements des pixels

Exercice 2.3 (Devine le résultat)
Que se passe-t-il lorsqu’on exécute les scripts suivants ?

Script 1 Script 2
clear all clear all

A=imread('lena512.bmp'); A=imread('lena.jpg');
colormap(gray(256)); [row,col]=size(A);

subplot(1,2,1) for j=1:col

imshow(uint8(A)); A=[A(:,2:col) A(:,1)];
title ( "Originale" ); imshow(A,gray(256));
subplot(1,2,2) pause(0.001);
B=A'; end
imshow(uint8(B));
title ( "Transposee" );
imwrite(uint8(B),'exotransposee.jpg','jpg');

Correction

1. Dans le premier script la matrice B est la transposée de la matrice A : l’image obtenue est la symétrique par rapport à la
diagonale principale :

Originale Transposée
2. À chaque étape on enlève la première colonne et on la concatène comme dernière colonne. Le résultat donne un “film”
dans lequel l’image “sort” à gauche et “rentre” à droite.

Exercice 2.4 (Flip, H-V-stack, Zoom)

En utilisant une manipulation élémentaire de la matrice (sans faire de boucles et sans utiliser de fonctions prédéfinies)
obtenir les images de la figure 2.2. On pourra se baser sur le canevas suivant :
clear all
A=imread('lena512.bmp');
B= .... % construire B
subplot(1,2,1)
imshow(uint8(A));
title ( "Originale" );
subplot(1,2,2)
imshow(uint8(B));
title ( "Transformee" );

(a) Original (b) Flip vertical (c) Flip horizontale (d) Flip double (e) Hstack (f) Vstack (g) Zoom

F IGURE 2.2 – Manipulations élémentaires

© 2022-2023 G. Faccanoni 39
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

Correction
1. Flip vertical B = A(end:-1:1,:);
Dans Octave, cela correspond à B = flipud(A); (up to down), dans Matlab à B = flip(A);

2. Flip horizontale B = A(:,end:-1:1);

Dans Octave, cela correspond à B = fliplr(A); (left to right), dans Matlab à B = flip(A,2);

3. Flip double B = A(end:-1:1,end:-1:1);

4. Hstack B = [ A(:,end/2:end), A(:,end:-1:end/2)];

5. Vstack B = [ A(end/2:end,:); A(end:-1:end/2,:)];

6. Zoom B = A(300:400,250:350);

Exercice 2.5 (Bandes)

On veut représenter un ensemble de 4 raies verticales blanches et 4 raies verticales noires alternées sur une image M × N
(M = N = 256).
1. Quelle est la taille en pixels de chaque raie blanche et chaque raie noire ?
2. Une première idée d’algorithme pour construire une telle image consiste d’abord a construire une raie blanche,
puis une raie noire puis former l’ensemble de l’image par juxtaposition.
³ ³ ´´
v( j −u)
3. L’image définie par b i , j = 256
2 1 + cos 2π N est aussi une image de raies. Comment choisir v et u pour que
l’image ainsi formée coïncide à peu prés avec l’image recherchée ?
4. On souhaite malgré tout une image binaire, à savoir une image pour laquelle les pixels n’ont que³ deux valeurs
´
v( j −u)
possibles : 0 ou 255. Modifier l’algorithme précédent en mettant en blanc les pixels pour lesquels cos 2π N ≥0
³ ´
v( j −u)
et en noir ceux pour lesquels cos 2π N < 0.
5. Décaler l’image vers la droite d’une demi-raie (i.e. à gauche de l’image et tout à droite il y a une demi-raie noire).

Correction
256 28
1. Chaque raie occupe 256 × 32 pixels. En effet, chaque raie occupe 256 lignes pour 8 = 23
= 25 = 32 colonnes.

2. clear all; clc;

N = zeros(256,32);
B = 255-N;
A = [B,N,B,N,B,N,B,N];
imshow(uint8(A));
%imwrite(uint8(A),'exoRaies1.jpg','jpg');

3. Il faut que b 16, j = 0, b 48, j = 1, etc.

clear all; clc;

B = 255*ones(256);
v = 32/8;
u = 32/2;
for j = 1:256
B(:,j) = 1+cos(2*pi*v*(j-u)/256);
end
% on se ramene a [0;255]
m = min(B(:));
M = max(B(:));
B = 255/(M-m).*(B-m);
imshow(uint8(B));
imwrite(uint8(B),'exoRaies2.jpg','jpg');

4.

40 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.2 Manipulations élémentaires : déplacements des pixels

clear all; clc;

B = 255*ones(256);
v = 32/8;
u = 32/2;
for j = 1:256
B(:,j) = cos(2*pi*v*(j-u)/256)> = 0;
end
% on se ramene a [0;255]
B = 255*B;
imshow(uint8(B));
imwrite(uint8(B),'exoRaies3.jpg','jpg');

5. clear all; clc;

N = zeros(256,32);
B = 255-N;
A = [B,N,B,N,B,N,B,N];
B = [A(:,16:end),A(:,1:15)];
imshow(uint8(B));
imwrite(uint8(B),'exoRaies4.jpg','jpg');

Exercice 2.6 (Réduction de la résolution)

Une matrice de taille 29 × 29 contient 218 entiers, ce qui prend pas mal de place en mémoire. On s’intéresse à des
méthodes qui permettent d’être plus économique sans pour cela diminuer la qualité esthétique de l’image. Afin de
réduire la place de stockage d’une image, on peut réduire sa résolution, c’est-à-dire diminuer le nombre de pixels. La
façon la plus simple d’effectuer cette réduction consiste à supprimer des lignes et des colonnes dans l’image de départ.
Les figures suivantes montrent ce que l’on obtient si l’on retient une ligne sur 2k et une colonne sur 2k ce qui donne une
matrice 29−k × 29−k . Appliquer cette transformation pour obtenir l’une des images suivantes :

(a) Originale (k = 1) (b) k = 2 (c) k = 3 (d) k = (e) k =

4 5

F IGURE 2.3 – Résolution

Correction subplot(1,8,k)
clear all E = A(1:2^k:row,1:2^k:col);
imshow(uint8(E));
A = double(imread('lena512.bmp')); title(['k = ' num2str(k)]);
colormap(gray(256)); file = strcat('exo2E', num2str(k) ,'.jpg')
%imwrite(uint8(E),file,'jpg');
[row,col] = size(A) end

for k = 1:8

© 2022-2023 G. Faccanoni 41
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

2.2.1 Transformation du photomaton

On part d’un tableau n × n, avec n pair, chaque élément du tableau représente un pixel. À partir de cette image on calcule
une nouvelle image en déplaçant chaque pixel selon une transformation, appelée transformation du photomaton.
On découpe l’image de départ selon des petits carrés de taille 2 × 2. Chaque petit carré est donc composé de quatre pixels. On
envoie chacun de ces pixels à quatre endroits différents de la nouvelle image : le pixel en haut à gauche reste dans une zone
en haut à gauche, le pixel en haut à droite du petit carré, est envoyé dans une zone en haut à droite de la nouvelle image,. . .

j =1 j =2 j ′ =1 j ′ =2

i =1 A B a b i ′ =1 A a B b

i =2 C D c d i ′ =2

C c D d

Par exemple le pixel en position (1, 1) (symbolisé par la lettre D) est envoyé en position (4, 4).
Explicitons ce principe par des formules. Pour chaque couple (i , j ), on calcule son image (i ′ , j ′ ) par la transformation du
photomaton selon les formules suivantes :
k+1
⋆ Si l’indice k est impair alors k ′ = 2 .
k+n
⋆ Si l’indice k est pair alors k ′ = 2 .
Voici un exemple d’un tableau 4 × 4 avant (à gauche) et après (à droite) la transformation du photomaton.

1 2 3 4 1 3 2 4
5 6 7 8 9 11 10 12
9 10 11 12 5 7 6 8
13 14 15 16 13 15 14 16

Exercice 2.7 (Photomaton)

1. Écrire une fonction photomaton(tableau) qui renvoie le tableau calculé après transformation. Par exemple le
tableau de gauche est transformé en le tableau de droite.

1 2 3 4 1 3 2 4
5 6 7 8 9 11 10 12
9 10 11 12 5 7 6 8
13 14 15 16 13 15 14 16

2. Écrire une fonction photomaton_iterer(tableau,k) qui renvoie le tableau calculé après k itérations de la
transformation du photomaton.
3. Expérimenter pour différentes valeurs de la taille n, afin de voir au bout de combien d’itérations on retrouve
l’image de départ.

Correction
Dans le fichier photomaton.m on écrit la fonction photomaton qui prend une matrice carrée de dimension n pair et renvoi
une matrice de même dimension qui a subi la transformation. Pour le calcul des nouveaux indices on écrit une sous-fonction
photomatonIndice :

function B = photomaton(A,n) wi = photomatIndice(i,n);

B = zeros(n); wj = photomatIndice(j,n);
for i = 1:n B(wi,wj) = A(i,j);
for j = 1:n end

42 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.2 Manipulations élémentaires : déplacements des pixels

end /2)
end w = floor((k+1)/2);
else
w = floor((k+d)/2);
function w = photomatIndice(k,d) end
if rem(k+1,2) = = 0 %(k+1)/2 = = floor((k+1) end

On teste la fonction d’abord sur la matrice 4 × 4 donnée puis sur une image :

’Script photomatonTEST.m’
clear all; clc;

% % TEST
% A = [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]
% B = photomaton(A,4)

A = imread('lena512.bmp');
[r,c] = size(A);
colormap(gray(256));

subplot(1,2,1)
imshow(uint8(A));
title ( "Original" );

subplot(1,2,2)
B = photomaton(A,r);
imshow(uint8(B));
title ( "Photomaton" );
%imwrite(uint8(B),'exoPhot1.jpg','jpg');

On boucle :

’Script photomatonTEST2.m’
clear all; clc;

A = imread('section1-original.png');
[r,c] = size(A);
colormap(gray(256));

subplot(3,3,1)
imshow(uint8(A));
for k = 2:9
subplot(3,3,k)
B = photomaton(A,r);
imshow(uint8(B));
A = B;
end

© 2022-2023 G. Faccanoni 43
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

2.2.2 Transformation du boulanger

Cette transformation s’appelle ainsi car elle s’apparente au travail du boulanger qui réalise une pâte feuilletée. L’image est
tout d’abord aplatie, coupée en deux puis la partie droite est placée sous la partie gauche en la faisant tourner de 180°.
On part d’un tableau n × n, avec n pair dont chaque élément représente un pixel. On va appliquer deux transformations
élémentaires à chaque fois :
⋆ Étirer : les deux premières lignes (chacune de longueur n) produisent une seule ligne de longueur 2n en mixant les
valeurs de chaque ligne en alternant un élément du haut, un élément du bas.

étirer

Voici comment deux lignes se mixent en une seule :

A B
A a B b
a b

Un élément en position (i , j ) du tableau de départ correspond à un élément (i /2, j /2) (si j est pair) ou bien ((i +1)/2, j /2)
(si j est impair) du tableau d’arrivée.
Exemple. Voici un tableau 4 × 4 à gauche, et le tableau étiré 2 × 8 à droite. Les lignes 0 et 1 à gauche donnent la ligne 0 à
droite. Les lignes 2 et 3 à gauche donne la ligne 1 à droite.

1 2 3 4
5 6 7 8 1 5 2 6 3 7 4 8
9 10 11 12 9 13 10 14 11 15 12 16
13 14 15 16

n
⋆ Replier : la partie droite d’un tableau étiré est retournée, puis ajoutée sous la partie gauche. Partant d’un tableau 2 × 2n
on obtient un tableau n × n.

A G
replier
A G
Pour 0 ≤ i < n2 et 0 ≤ j < n les éléments en position (i , j ) du¡tableau sont conservés. Pour n2 ≤ i < n et 0 ≤ j < n un
n
¢
élément du tableau d’arrivée (i , j ), correspond à un élément 2 − i − 1, 2n − 1 − j du tableau de départ.
Exemple. À partir du tableau étiré 2 × 8 à gauche, on obtient un tableau replié 4 × 4 à droite.

1 5 2 6
1 5 2 6 3 7 4 8 9 13 10 14
9 13 10 14 11 15 12 16 16 12 15 11
8 4 7 3

La transformation du boulanger est la succession d’un étirement et d’un repliement. Partant d’un tableau n × n on obtient
encore un tableau n × n.

Exercice 2.8 (Boulanger)

1. Écrire une fonction boulanger(tableau) qui renvoie le tableau calculé après transformation. Par exemple le
tableau de gauche est transformé en le tableau de droite.

1 2 3 4 1 5 2 6
5 6 7 8 9 13 10 14
9 10 11 12 16 12 15 11
13 14 15 16 8 4 7 3

2. Écrire une fonction boulanger_iterer(tableau,k) qui renvoie le tableau calculé après k itérations de la
transformation du photomaton.

44 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.2 Manipulations élémentaires : déplacements des pixels

3. Vérifier que pour n = 256 au bout de 17 itérations on retrouve l’image de départ.

Correction
Dans le fichier boulanger.m on écrit la fonction boulanger qui prend une matrice carrée de dimension n pair et renvoi
une matrice de même dimension qui a subi la transformation.

function C=boulanger(A,n) B(fix((i+1)/2),2*j) = A(i,j);

B = zeros(fix(n/2),2*n); end
% etirer end
for i = 1:n end
for j = 1:n % couper
if fix((i+1)/2)==(i+1)/2 C = [B(:,1:n); B(end:-1:1,end:-1:n+1)];
B(fix((i+1)/2),2*j-1) = A(i,j); end
else

On teste nos fonctions d’abord sur la matrice 4 × 4 donnée puis sur une image :

’Script boulangerTEST.m’
clear all; clc;

% TEST
A = [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]
B = boulanger(A,4)

A = imread('lena512.bmp');
[r,c] = size(A);
colormap(gray(256));

subplot(1,2,1)
imshow(uint8(A));
title ( "Original" );
subplot(1,2,2)
B=boulanger(A,r);
imshow(uint8(B));
title ( "Boulanger" );
%imwrite(uint8(B),'exoBou1.jpg','jpg');

On boucle :

’Script boulangerTEST2.m’
clear all; clc;

A = imread('section1-original.png');
[r,c] = size(A);
colormap(gray(256));

subplot(4,5,1)
imshow(uint8(A));
for k = 2:18
disp(k)
subplot(4,5,k)
B = boulanger(A,r);
imshow(uint8(B));
A = B;
end

© 2022-2023 G. Faccanoni 45
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

2.3 Manipulations élémentaires : modification des valeurs des pixels

On s’intéresse aux traitements ponctuels des images numériques qui consistent à faire subir à chaque pixel une correction ne
dépendant que de sa valeur.

Exercice 2.9 (Censure, Bords)

En utilisant une manipulation élémentaire de la matrice (sans faire de boucles et sans utiliser de fonctions prédéfinies)
obtenir les images de la figure 2.4 (pour la dernière image, regarder la doc de meshgrid sinon utiliser deux boucles). On
pourra se baser sur le canevas suivant :
clear all
A=imread('lena512.bmp');
B= .... % construire B
subplot(1,2,1)
imshow(uint8(A));
title ( "Originale" );
subplot(1,2,2)
imshow(uint8(B));
title ( "Transformee" );

(a) Original (b) Effacer (c) Bord carré (d) Bord cercle

F IGURE 2.4 – Manipulations élémentaires

Correction
1. Effacer
B = A; B(250:300,220:375)= 0;

2. Bord carré
B = A; [r,c] = size(A); B([1:20,r-20:r],:)= 0; B(:,[1:20,c-20:c])= 0;

3. Bord cercle
Avec meshgrid (sans boucles) :
[l,c] = size(A); r = min(l,c);
[xx,yy] = meshgrid(1:l,1:c);
M = (xx-l/2).^2+(yy-c/2).^2 >= (r/2).^2;
B = A; B(M) = 255;

Avec boucles :
[r,c]=size(A);
B=255*ones(r,c);
for i=1:r
for j=1:c
if (i-r/2)^2+(j-c/2)^2 <= (r/2)^2
B(i,j)=A(i,j);
end
end
end

Pour comprendre ce que renvoie la fonction meshgrid,

considérons les points suivants positionnés sur une grille 2
cartésienne : 1

3 4 5

46 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.3 Manipulations élémentaires : modification des valeurs des pixels

Il s’agit de l’ensemble de points x = 3:5 ;

y = 1:2 ;
{ (3, 1); (4, 1); (5, 1); (3, 2); (4, 2); (5, 2) } . [xx,yy] = meshgrid(x,y)

Pour générer deux matrices qui contiennent respective- xx =

ment les abscisses et les ordonnées correspondantes, on 3 4 5
utilise la fonction meshgrid. 3 4 5

yy =
1 1 1
2 2 2

Exercice 2.10 (Statistiques, Histogramme)

On considère toutes les valeurs des pixels de l’image. Le calcule de la moyenne, de la variance, de l’écart type et de la
médiane ainsi que la visualisation de l’histogramme donnent de nombreuses informations sur l’image. L’histogramme
h d’une image associe à chaque valeur d’intensité de gris le nombre de pixels prenant cette valeur. La détermination de
l’histogramme est donc réalisée en comptant le nombre de pixel pour chaque intensité de l’image. Un histogramme
indique, pour chaque valeur entre le noir (0) et le blanc (255), combien il y a de pixels de cette valeur dans l’image ; en
abscisse (axe x) le niveau de gris (de 0 à 255) ; en ordonnée (axe y) le nombre de pixels. Les pixels sombres apparaissent
à gauche de l’histogramme, les pixels clairs à droite de l’histogramme et les pixels gris au centre de l’histogramme.
1. Calculer la moyenne, l’écart type et la médiane d’une image.
2. Tracer son histogramme.

Correction

1. clear all
clc
A=double(imread('lena512.bmp')); % pour travailler avec des valeurs de type double
moyenne=mean(A(:))
mediane=median(A(:))
ecart_type=std(A(:))

2. clear all
A=imread('lena512.bmp');
subplot(2,1,1)
colormap(gray(256));
A=double(A);
imshow(uint8(A));
subplot(2,1,2)
x=[0:255];
h=[];
for i = x
h=[h, sum(A(:)==i)];
end
moy=mean(A(:));
med=median(A(:));
plot(x,h,'-',[moy moy],[0,max(h)],[med med
],[0,max(h)])
title([ "Moy=",num2str(moy), " - Mediane=",
num2str(med) ])
axis([0 255 0 max(h)],"square")
grid()

print('lenahist.png', '-dpng');

Remarque : avec Octave, la fonction prédéfinie imhist calcule l’histogramme de l’image mais elle nécessite l’installation du
package image. C’est pourquoi on a préféré ici construire l’histogramme manuellement :

© 2022-2023 G. Faccanoni 47
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

clear all stem(nb,counts);

figure(1,"position",get(0,"screensize")) axis([0 255 0 max(counts)])
A=imread('lena512.bmp'); grid()

pkg load image % Rque : h'=counts

[counts, nb]=imhist(A,gray(256));

⋆
Une image trop lumineuse aura une moyenne et une médiane supérieures à 200 ainsi qu’un histogramme déplacés
vers la droite (on dit que l’image est brûlée). Une image peu lumineuse aura une moyenne et une médiane inférieures
à 100 ainsi qu’un histogramme déplacés vers la gauche (on dit que l’image est bouchée).
⋆ Une image à faible contraste aura un petit écart type ainsi qu’un histogramme concentré autour de la médiane. Une
image à fort contraste aura un grand écart type ainsi qu’un histogramme étalé.
Suivant les informations données par la moyenne, l’écart type et l’histogramme, on pourra améliorer visuellement l’image en
modifiant la valeur de chaque pixel par une fonction f qu’on appliquera à tous les pixels. Pour cela, on pourra se baser sur le
canevas suivant :
clear all
A=double(imread("lena512.bmp")); % utiliser double pour avoir des calculs precis

f=@(x) .... ; % fonction vectorisee a completer

B=f(A);
subplot(1,3,1)
plot([0:255], f([0:255]), "b-",[0:255], [0:255], "r--"); % f et identite
axis([0 255 0 255],"square");
title("f vs identite")
subplot(1,3,2)
imshow(uint8(A));
title ( "Originale" );
subplot(1,3,3)
imshow(uint8(B));
title ( "Transformee" );

y
Par exemple, la fonction suivante, appliquée à chaque pixel 255
d’une image, éclaircira l’image :

f : [0; 255] → [0; 255]

a i j + 255
a i j 7→
2 0
0 255 a i j
Dans le canevas, il suffit de définir
f = @(g) (g+255)/2;

et on obtient

Exercice 2.11 (Négatif )

Pour obtenir l’image en négatif de Léna il suffit de prendre le complémentaire par rapport à 255

f : [0; 255] → [0; 255]

x 7→ 255 − x

48 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.3 Manipulations élémentaires : modification des valeurs des pixels

Appliquer cette transformation pour obtenir l’image 2.5b.

(a) Originale (b) Négatif

F IGURE 2.5 – Négatif

Correction
f = @(x)255-x;

Exercice 2.12 (Luminosité)

1. Pour augmenter la luminosité, il suffit d’ajouter une valeur fixe d > 0 à tous les niveaux. Pour diminuer la
luminosité il faudra au contraire soustraire une valeur fixe d > 0 à tous les niveaux.
Appliquer la transformation ci-dessous pour augmenter la luminosité en ajoutant la valeur fixe d = 50 à tous les
niveaux de gris et obtenir l’image 2.6b.
y
255
f : [0; 255] → [0; 255]
(
x +d si x ≤ 255 − d , d
x 7→
255 sinon.
0
0 255 − d 255 x
2. Il est très mauvais d’augmenter ainsi la luminosité : avec un décalage de d , il n’existera plus aucun point entre 0
et d et les points ayant une valeur supérieure à 255 − d deviendront des points parfaitement blancs, puisque la
valeur maximale possible est 255. La nouvelle image contient des zones dites “brûlées”.
y
Plutôt que d’utiliser la fonction donnée, il vaut mieux 255
utiliser une fonction bijective de forte croissance au voi-
sinage de 0 et de très faible croissance au voisinage de
255, comme sur le graphe ci-contre.
Appliquer cette transformation pour obtenir l’image 2.6c.
0
0 255 x

(a) Originale (b) Première méthode (c) Deuxième méthode

F IGURE 2.6 – Modification de la luminosité

Correction
1. f = @(g)min(255,g+50);

2. On cherche une fonction f continue telle que f (0) = 0, f (255) = 255 et f (x) > x si x ∈]0; 255[.
p
⋆ On peut considérer par exemple la fonction f (x) = 255x. La fonction f s’écrit f = @(x)sqrt(255*x);

© 2022-2023 G. Faccanoni 49
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

⋆ On peut considérer le quart de cercle de centre (255, 0) et rayon 255. Il a pour équation
p (x − 255)2 + (y − 0)2 = 2552 .
Le quart de cercle qui nous intéresse est la partie qui a pour équation f (x) = 255 − (x − 255)2 . La fonction f
2

s’écrit f = @(x)sqrt(255^2-(x-255).^2);
⋆ On peut chercher cette fonction sous la forme d’une parabole, donc d’équation y = ax 2 + bx + c. L’inégalité se
traduit par a < 0. Les deux conditions f (0) = 0 et f (255) = 255 ne donnent que deux équations linéaires en les
trois inconnues a, b, c, il faut donc ajouter une autre condition pour fixer une parabole. On peut par exemple
b
choisir la parabole dont le sommet se trouve en x = 255, i.e. la condition − 2a = 255. On a alors les trois conditions

c =0
    
 0 0 1 a 0
2552 a + 255b + c = 255 ⇐⇒  2552 255 1 b  = 255
2 × 255 1 0 c 0

2 × 255a + b = 0


1
On peut résoudre le système à la main et on trouve c = 0, a = − 255 et b = 2 ou demander à Matlab/Octave : [0 0
1;255^2 255 1;2*255 1 0]\[0;255;0] et finalement la fonction f s’écrit f = @(x)x.*(255*2-x)/255;

Exercice 2.13 (Contraste)

1. On se donne une image présentant un histogramme concentré dans l’intervalle [m; M ]. Les valeurs m = min(A) et
M = max(A) correspondent aux niveaux de gris extrêmes présents dans cette image. Le recadrage de dynamique
consiste à étendre la dynamique de l’image transformée à l’étendue totale [0; 255] (Histogram Stretching). La
transformation de recadrage est donc une application affine qui s’écrit :

y
f : [0; 255] → [0; 255]
255

0
 si x ≤ m
255−0
x 7→ M −m (x − m) + 0 si m < x < M

255 si x ≥ M


0
Il s’agit de la droite qui passe par les points (m, 0) et 0 m M 255 x
(M , 255) sur [m, M ].
Appliquer cette transformation pour obtenir l’image 2.7b.
2. On peut augmenter le contraste (Contrast Stretching) en “mappant” par une fonction linéaire par morceaux.
y
255
f : [0; 255] → [0; 255]

0
 si x ≤ a
255−0
x 7→ b−a (x − a) + 0 si a < x < b

255 si x ≥ b

0
0 a b 255 x
Appliquer cette transformation avec a = 64 et b = 192 pour obtenir l’image 2.7c.
Un cas particulier s’obtient lorsque a = b (Contrast Thresholding ou seuillage). Appliquer cette transformation
pour obtenir l’image 2.7d
Avec cette transformation, les points les plus blancs au- y
ront une valeur égale à b et il n’existera plus aucun point 255
entre b et 255. De même, les points ayant une valeur com-
prise entre 0 et a deviendront noirs, puisque la valeur m
minimale est 0. Il y aura donc là perte d’informations.
Ainsi il est plus judicieux d’adoucir la courbe comme dans
0
l’image ci-contre : 0 m 255 x
3. On peut rehausser le contraste en “mappant” par une fonction linéaire par morceaux

50 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.3 Manipulations élémentaires : modification des valeurs des pixels

y y
255 255
f : [0; 255] → [0; 255] b>a
b

 x si x ≤ a


a b<a
x 7→
 (255 − b)x + 255(b − a)
si x ≥ a


255 − a 0 0
0 a 255 x 0 a 255 x
Appliquer cette transformation avec a = 100 et b = 200 (dilatation de la dynamique des zones claires) et comparer
avec a = 100 et b = 50 (dilatation de la dynamique des zones sombres). Appliquer ces transformations pour obtenir
les images 2.7e et 2.7f.
4. On pourrait vouloir mettre en avant certains niveaux de gris dans un certain intervalle. Cette approche peut
prendre deux formes : le gray level slicing without background et le gray level slicing with background. Dans la
première approche, une grande valeur est donnée aux pixels dont le niveau de gris appartient à la plage choisie et
les autres sont remplacés par 0. Dans la deuxième approche on ne modifie que les pixels à mettre en valeur. Cela
correspond aux deux fonctions suivantes :
y
255
f : [0; 255] → [0; 255]
(
0 si x ≤ a ou x ≥ b
x 7→
255 si a ≤ x ≤ b
0
0 a b 255 x
y
255
f : [0; 255] → [0; 255]
(
x si x ≤ a ou x ≥ b
x 7→
255 si a ≤ x ≤ b
0
0 a b 255 x
Appliquer ces deux transformations avec a = 100 et b = 155 pour obtenir les images 2.7g et 2.7h .
5. On peut modifier le contraste en “mappant” par une fonction croissante plus rapidement ou plus lentement que
la fonction identité i : [0; 255] → [0; 255], i (x) = x. Par exemple, la fonction suivante modifie la caractéristique de
gamma (plus clair si 0 < a < 1, plus foncé si a > 1).
y
255

f : [0; 255] → [0; 255]

³ x ´a
x 7→ 255
255
0
0 255 x
Appliquer cette transformation avec a = 3 et a = 1/3 pour obtenir les images 2.7i et 2.7j.

© 2022-2023 G. Faccanoni 51
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

(a) Originale (b) Histogram Stret- (c) Contrast Stretching (d) Contrast Threshol- (e) Dilatation de la dy-
ching ding namique des zones
claires

(f) Dilatation de la dy- (g) Gray level slicing (h) Gray level slicing (i) Non linear (j) Non linear
namique des zones without back- with background
sombres ground

F IGURE 2.7 – Modification du contraste

Correction
1. Histogram Stretching :
Tout d’abord il faut calculer m et M . On écrit 4
soit M=max(A(:)); m=min(A(:));
soit M=max(max(A)); m=min(min(A));
Ensuite il faut définir la fonction f :
f = @(x)(x<=m)*0 + (x>=M)*255 + (x>m).*(x<M).*( 255/(M-m)*(x-m));
Cependant, on peut juste définir f comme
f = @(x)255*(x-m)/(M-m);
En effet, on obtiendra f (x) < 0 si x < m et f (g ) > 255 si x > 255 et comme on applique la fonction uint8 à B , les valeurs
négatifs seront considérés égales à 0 et les valeurs supérieurs à 255 seront considérés égales à 255.

2. Contrast Stretching :
a=64; b=192; f = @(x)(x<=a)*0 + (x>=b)*255 + (x>a).*(x<b).*( 255/(b-a)*(x-a));
ou, pour la même raison que le point précédent,
a=64; b=192; f = @(x)255/(b-a)*(x-a);
Si a = m et b = M on retrouve l’Histogram Stretching.

Contrast Thresholding :
a=128; f = @(x)(x<=a)*0 + (x>=a)*255;

Amélioration :
On cherche une fonction f continue telle que f (0) = 0, f (m) = m, f (255) = 255, f (x) < x si x < m et f (x) > x si x > m.
⋆ On peut considérer les deux quarts de cercle suivants. Le premier est le cercle de centre (0, m) et rayon
p m qui a
pour équation (x − 0)2 + (y − m)2 = m 2 et la partie qui nous intéresse a pour équation f (x) = m − m 2 − x 2 . Le
deuxième est le cercle de centre (255, m) et rayon 255 − mpqui a pour équation (x − 255)2 + (y − m)2 = (255 − m)2 et
la partie qui nous intéresse a pour équation f (x) = m + (255 − m)2 − (x − 255)2 .
La fonction f s’écrit f = @(x)cb(x).*(x<m)+ch(x).*(x>=m); avec cb = @(x)m-sqrt(m^2-x.^2); et
ch = @(x)m+sqrt((255-m)^2-(x-m).^2);

4. Pour comprendre voir cette exemple :

A=[1 2;3 4]
A(:)
min(A)

52 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.3 Manipulations élémentaires : modification des valeurs des pixels

y
255

0
0 m 255 x

⋆ On peut choisir de définir f par morceaux avec deux paraboles :

(
p 1 (x) si x ≤ m,
f (x) =
p 2 (x) si x ≥ m.

On a p 1 (0) = 0, p 1 (m) = m, p 1 est convexe mais il faut ajouter une condition pour fixer cette parabole, par exemple
on demande à ce que le sommet se trouve en 0. On obtient ainsi p 1 (x) = x 2 /m :
p1=@(x)(x.^2/m);
On a p 2 (255) = 255, p 2 (m) = m, p 2 est concave mais il faut ajouter une condition pour fixer cette parabole, par
exemple on demande à ce que le sommet se trouve en 255. Si on écrit la parabole sous la forme p 2 (x) = ax 2 +bx +c,
on doit résoudre le système linéaire

2

 255 a + 255b + c = 255 2552
    
 255 1 a 255
m 2 a + mb + c = m ⇐⇒  m2 m 1  b  =  m 
2 × 255 1 0 c 0

2 × 255a + b = 0


On peut bien sûr résoudre ce système à la main, mais on peut demander à Matlab/Octave de le faire pour nous :
sol=[255^2 255 1; m^2 m 1; 2*255 1 0]\[255;m;0];
a=sol(1);
b=sol(2);
c=sol(3);
p2=@(x)(a*x.^2+b*x+c);
et enfin on a la fonction
f=@(x)p1(x).*(x<=m)+ p2(x).*(x>m)

3. Dilatation de la dynamique des zones claires :

a=100; b=200; f = @(x)(x<=a)*(b/a).*x + (x>a).*( ((255-b)*x+255*(b-a))/(255-a));

Dilatation de la dynamique des zones sombres :

a=100; b=50; f = @(x)(x<=a)*(b/a).*x + (x>a).*( ((255-b)*x+255*(b-a))/(255-a));

4. Gray level slicing without background :

a=85; b=170; f = @(x)(x<=a)*0 + (x>=b)*0 + (x>a).*(x<b)*255 ;
ou, plus simplement
a=85; b=170; f = @(x)(x>a).*(x<b)*255 ;

Gray level slicing with background :

a=85; b=170; f = @(x)(x<=a).*x + (x>=b).*x + (x>a).*(x<b)*255 ;
Attention à ne pas écrire f = @(x)(x<=a).*(x>=b).*x + (x>a).*(x<b)*255 ; car le produit (x<=a).*(x>=b)
correspond à x ≤ a ET x ≥ b alors que nous voulons x ≤ a OU x ≥ b.

5. Non linear :
a=3; f = @(x)255*(x/255).^a;
a=1/3; f = @(x)255*(x/255).^a;

Exercice 2.14 (Quantification)

Une autre façon de réduire la place mémoire nécessaire pour le stockage consiste à utiliser moins de nombres entiers
pour chaque valeur. On peut par exemple utiliser uniquement des nombres entiers entre 0 et 3, ce qui donnera une
image avec uniquement 4 niveaux de gris. Une telle opération se nomme quantification.
On peut effectuer une conversion de l’image d’origine vers une image avec 28−k niveaux de valeurs en effectuant les
remplacements suivant : tous les valeurs entre 0 et 2k sont remplacées par la valeur 0, puis tous les valeurs entre 2k et

© 2022-2023 G. Faccanoni 53
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

2k+1 sont remplacées par la valeur 2k etc. Appliquer cette transformation pour obtenir les images suivantes :

(a) 16 niveaux de gris (k = 4) (b) 8 niveaux de gris (k = 5) (c) 4 niveaux de gris (k = 6) (d) 2 niveaux de gris (k = 7)

F IGURE 2.8 – Quantification

Correction subplot(2,2,3)
clear all hist(A(:),0:255);
subplot(2,2,2)
A = imread('lena512.bmp'); B = fix(A/2^k)*2^k; %idivide(A,2^k)*2^k;
colormap(gray(256)); imshow(uint8(B));
A = double(A); title ( strcat("Quantifiee, k = ",num2str(k)) )
[row,col] = size(A) ;
subplot(2,2,4)
for k = [4,5,6,7] hist(B(:),0:255);
figure(k) %imwrite(uint8(B),strcat("exo3",num2str(k-3),".
subplot(2,2,1) jpg"),'jpg');
imshow(uint8(A)); end
title ( "Original" );

54 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.4 Détection des bords

2.4 Détection des bords

Jusqu’à maintenant, le traitement de chaque pixel d’une image était indépendant des pixels voisins. Il s’agit maintenant
d’aborder des algorithmes dont le résultat est lié aux pixels voisins.
Afin de localiser des objets dans les images, il est nécessaire de détecter les bords de ces objets. Ces bords correspondent à
des zones de l’image où les valeurs des pixels changent rapidement. C’est le cas par exemple lorsque l’on passe du chapeau
(qui est clair, donc avec des valeurs grandes) à l’arrière plan (qui est sombre, donc avec des valeurs petites). Afin de savoir
si un pixel avec une valeur est le long d’un bord d’un objet, on prend en compte les valeurs de ses quatre voisins (deux
horizontalement et deux verticalement) en calculant une valeur b i , j suivant la formule
q
bi , j = (a i +1, j − a i −1, j )2 + (a i , j +1 − a i , j −1 )2

On peut remarquer que si b i , j = 0, alors on a a i −1, j = a i +1, j et a i , j −1 = a i , j +1 . Au contraire, si b i , j est grand, ceci signifie que
les pixels voisins ont des valeurs très différentes, le pixel considéré est donc probablement sur le bord d’un objet. Notons que
l’on utilise ici seulement 4 voisins, ce qui est différent du calcul de moyennes où l’on utilisait 8 voisins. Ceci est important
afin de détecter aussi précisément que possible les bords des objets.

Exercice 2.15 (Détection des bords)

Appliquer cette transformation pour obtenir les images 2.9b et 2.9c. Pour améliorer le rendu, ramener les niveaux de
gris à l’intervalle [0; 255] par une transformation affine et obtenir l’image 2.9d.

(a) Original (b) Contour (c) Contour Négatif (d) Contour Négatif
(normalisé)

F IGURE 2.9 – Masques

Correction
clear all

A = double(imread('lena512.bmp'));
colormap(gray(256));

[row,col] = size(A);
B = zeros(row,col);

I = 2:row-1;
J = 2:col-1;
B(I,J) = sqrt( (A(I+1,J)-A(I-1,J)).^2+(A(I,J+1)-A(I,J-1)).^2 );

% equivaut aux boucles suivantes :

%for i = 2:row-1
% for j = 2:col-1
% B(i,j) = sqrt( (A(i+1,j)-A(i-1,j))^2+(A(i,j+1)-A(i,j-1))^2 );
% end
%end

subplot(1,4,1)
imshow(uint8(A));
title ( "Original" );
%imwrite(uint8(A),'exo6orig.jpg','jpg');
subplot(1,4,2)
imshow(uint8(B));
title ( "Contour" );
%imwrite(uint8(B),'exo6grad.jpg','jpg');
subplot(1,4,3)

© 2022-2023 G. Faccanoni 55
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

imshow(uint8(255-B));
title ( "Negatif du contour" );
%imwrite(uint8(255-B),'exo6gradNeg.jpg','jpg');
subplot(1,4,4)
% on se ramene a [0;255]
m = min(min(B(I,J)))
M = max(max(B(I,J)))
B = 255/(M-m).*(B-m);
imshow(uint8(255-B));
title ( "Negatif du contour (normalise)" );
%imwrite(uint8(255-B),'.jpg','jpg');

2.5 Masques (filtrage par convolution)

Exercice 2.16 (Enlever du bruit par moyennes locales)
Les images sont parfois de mauvaise qualité. Un exemple typique de défaut est le bruit qui apparaît quand une photo est
sous-exposée, c’est-à-dire qu’il n’y a pas assez de luminosité. Ce bruit se manifeste par de petites fluctuation aléatoires
des niveaux de gris. La figure 2.10b montre une image bruitée obtenue par modification de l’image 2.10a par les
instructions
sigma=50;
B=A+randn(size(A))*sigma;

Afin d’enlever le bruit dans les images, il convient de faire subir une modification aux valeurs de pixels. L’opération la
plus simple consiste à remplacer la valeur a i j de chaque pixel par la moyenne de a i j et des 8 valeurs a i −1, j −1 , a i −1, j ,
a i −1, j +1 , a i , j −1 , a i , j +1 , a i +1, j −1 , a i +1, j , a i +1, j +1 des 8 pixels voisins de a i j .
En effectuant cette opération pour chaque pixel, on supprime une partie du bruit, car ce bruit est constitué de fluc-
tuations aléatoires, qui sont diminuées par un calcul de moyennes. Appliquer cette transformation pour obtenir
l’image 2.10c Cependant, tout le bruit n’a pas été enlevé par cette opération. Afin d’enlever plus de bruit, on peut
moyenner plus de valeurs autour de chaque pixel. Le moyennage des pixels est très efficace pour enlever le bruit
dans les images, malheureusement il détruit également une grande partie de l’information de l’image. on peut en effet
s’apercevoir que les images obtenues par moyennage sont floues (c’est justement l’effet recherché lors de la pixellisation).
Ceci est en particulier visible près des contours, qui ne sont pas nets.
Afin de réduire ce flou, on peut remplacer la moyenne par la médiane. Appliquer cette transformation pour obtenir
l’image 2.10d

(a) Originale (b) Image bruitée (c) Image débruitée (d) Image débruitée
par moyenne par médiane

F IGURE 2.10 – Dé bruitage

Correction for i=2:r-2

clear all, close all for j=2:c-2
sousB=B(i-1:i+1,j-1:j+1);
A=double(imread('lena512.bmp')); C(i,j)=mean(sousB(:));
[r,c]=size(A); D(i,j)=median(sousB(:));
end
sigma=50; end
B=A+randn([r,c])*sigma;
C=B; subplot(1,4,1)
D=B; imshow(uint8(A),[])
title("Image originale");

56 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.5 Masques (filtrage par convolution)

imwrite(uint8(A),strcat("exoAA.jpg"),'jpg'); imshow(uint8(C),[])
title("Image debruitee par moyenne");
subplot(1,4,2) imwrite(uint8(C),strcat("exoCC.jpg"),'jpg');
imshow(uint8(B),[])
title("Image bruitee"); subplot(1,4,4)
imwrite(uint8(B),strcat("exoBB.jpg"),'jpg'); imshow(uint8(D),[])
title("Image debruitee par mediane");
subplot(1,4,3) imwrite(uint8(D),strcat("exoDD.jpg"),'jpg');

Exercice 2.17 (Pixellisation)

La pixellisation est souvent utilisée pour flouter une image, en partie ou dans sa totalité. La zone pixelisée fait apparaître
des blocs carrés contenant n × n pixels de couleur uniforme. On observe la même chose en agrandissant une image
jusqu’au niveau du pixel, sauf que là, chaque carré représente 1 pixel.
Lorsqu’on diminue la résolution, les images sont de plus en plus petites (dans l’exercice à chaque étape on divisait par
deux le nombre de lignes et de colonnes de la matrice). Ici nous voulons obtenir des images qui ont toujours la même
dimension mais avec cet effet de pixellisation. Pour arriver à ce résultat, le principe de l’algorithme est de diviser l’image
en blocs carrés de largeur n pixels et remplacer tous les pixels du bloc par leur moyenne. On parcourt cette image non
pas pixel par pixel mais par blocs de n × n pixels.
Appliquer cette transformation avec n = 16 pour obtenir l’image 2.11b

(a) Originale (b) n = 16

F IGURE 2.11 – Pixellisation

Correction %mask = A(i:i+n-1,j:j+n-1);

clear all %B(i:i+n-1,j:j+n-1) = mean(mask(:));
end
A = double(imread('lena512.bmp')); end
colormap(gray(256));
[r,c] = size(A); subplot(1,2,1)
n = 16; imshow(uint8(A));
B = A; title ( "Original" );
for i=1:n:r-n+1 subplot(1,2,2)
for j=1:n:c-n+1 imshow(uint8(B));
% octave title (strcat([ "Pixellisation n=" num2str(n)]) );
B(i:i+n-1,j:j+n-1) = mean(A(i:i+n-1,j:j+n-1) %imwrite(uint8(B),'pixellisation.jpg','jpg');
(:));
%%%%%%% matlab

Convolution
Le principe du filtrage est de modifier la valeur des pixels d’une image, généralement dans le but d’améliorer son aspect. En
pratique, il s’agit de créer une nouvelle image en se servant des valeurs des pixels de l’image d’origine.
Un filtre est une transformation mathématique (appelée produit de convolution) permettant de modifier la valeur d’un pixel
en fonction des valeurs des pixels avoisinants, affectées de coefficients. Le filtre est représenté par un tableau (une matrice),
caractérisé par ses dimensions et ses coefficients, dont le centre correspond au pixel concerné. La somme des coefficients
doit faire 1.
Voici un exemple. Soit A une matrice et considérons le pixel a i , j et les 9 pixels voisins. On peut construire une matrice B où le
pixel b i , j est une combinaison linéaire (i.e. une moyenne pondérée) de tous ces pixels : 5

b i , j = c 1 a i −1, j −1 + c 2 a i −1, j + c 3 a i −1, j +1 + c 4 a i , j −1 + c 5 a i , j + c 6 a i , j +1 + c 7 a i +1, j −1 + c 8 a i +1, j + c 9 a i +1, j +1

On peut visualiser cela comme une masque : on considère A i −1,...,i +1 la sous-matrice carrée d’ordre 3 et le masque M
j −1,..., j +1

5. La relation ci-dessus se généralise aux masques de convolution 5 × 5, 7 × 7, etc.

© 2022-2023 G. Faccanoni 57
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

suivantes :
a i −1, j −1 a i −1, j a i −1, j +1 c1 c2 c3
A i −1,...,i +1 = a i , j −1 ai , j a i , j +1 M = c4 c5 c6
j −1,..., j +1 a i +1, j −1 a i +1, j a i +1, j +1 c7 c8 c9
Rappelons que les valeurs des composantes des pixels sont des nombres entiers compris entre 0 et 255. Si les nouvelles
valeurs ne sont plus des entiers, il faudra les arrondir. Il peut même arriver que la nouvelle valeur ne soit plus comprise entre
0 et 255, il faudra alors choisir si les écrêter ou appliquer une transformation affine.
De plus, on voit que la relation ne peut s’appliquer sur les bords de l’image. Il faut donc prévoir un traitement spécial pour
ces rangées. La solution la plus simple, consiste à remplir ces rangées avec un niveau constant (par exemple noir). Une autre
solution est de répliquer sur ces rangées les rangées voisines. Dans ces deux cas, l’information sur les bords est perdue. On
peut aussi décider de ne pas modifier ces rangées, ce que nous allons faire. En tout cas, on évite de réduire la taille de l’image.

Exercice 2.18 (Devine le résultat) ³0 0 0´

Appliquer une fois le masque suivant à la matrice A = 010 :
000

1/9 1/9 1/9

1/9 1/9 1/9
1/9 1/9 1/9

Correction
Seul le pixel en position (2, 2) sera modifié car on a décidé de ne pas modifier les pixels de bord :

1 X 3 1 1
b 2,2 = a i j = a 2,2 =
9 i , j =1 9 9

ainsi  
0 0 0
B = 0 1/9 0
0 0 0

Exercice 2.19 (Lissage, floutage, accentuation)

Écrire une fonction qui prend en entrée A (l’image) et M le masque (matrice d’ordre 3) et renvoie B l’image modifiée.
Avant d’afficher A ou B appliquer une transformation affine pour que les valeurs soient comprises entre 0 et 255.
1. Lissage (filtre passe-bas). On veut lisser les endroits à fort gradient pour enlever du bruit ou flouter une image. Le
principe de la plupart des méthodes de débruitage est simple : il consiste à remplacer la valeur d’un pixel par la
valeur moyenne des pixels voisins. La réduction du bruit s’accompagne d’une réduction de la netteté. Appliquer
50 fois le masque suivant pour obtenir l’image 2.12b.

1/9 1/9 1/9

1/9 1/9 1/9
1/9 1/9 1/9

2. Floutage d’une partie. Appliquer 50 fois le masque précédent à une sous-matrice bien choisie pour obtenir
l’image 2.12c.
3. Accentuation (filtre passe-haut). Appliquer 1 fois le masque suivant pour obtenir l’image 2.12d.

0 −1/2 0
−1/2 3 −1/2
0 −1/2 0

4. Rehaussement. L’objectif de rehaussement de contours est de rendre plus clairs et nets les contours en réduisant
la transition de l’intensité. Donc un opérateur de rehaussement vise à remplacer le pixel central par la somme des
différences avec ses voisins. Dans un premier temps, appliquer simplement le masque suivant à l’image de départ
ce qui permet de détecter les contours. Puis retrancher ce résultat multiplié par un coefficient α = 0.1 à l’image de
départ pour obtenir l’image 2.12e. Cette dernière opération permet dans le cas d’un échelon moue (contour pas

58 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.5 Masques (filtrage par convolution)

très distingué) de rehausser le contour si α est positif.

1/8 1/8 1/8

1/8 −1 1/8
1/8 1/8 1/8

5. Contours. Ci dessous on a indiqué deux masques. Appliquer le premier masque à l’image de départ pour obtenir
une image qu’on appellera X. Appliquer le deuxième masque à l’imageqde départ pour obtenir une image qu’on
appellera Y. Afficher l’image associée à la matrice G définie par Gi j = X2i j + Y2i j pour obtenir l’image 2.12f.

1/8 0 −1/8 1/8 1/4 1/8

1/4 0 −1/4 0 0 0
1/8 0 −1/8 −1/8 −1/4 −1/8

Même exercice avec les deux masques :

0 −1 0 0 0 0
0 0 0 −1 0 1
0 1 0 0 0 0

(a) Original (b) Floutage (c) Floutage sous- (d) Filtre passe-haut (e) Rehaussement (f) Contours
matrice

F IGURE 2.12 – Masques

Correction
D’abord on définit la fonction :
function Bnew=myflou(A,K,n)
[r,c]=size(A);
Bold=A;
Bnew=Bold;
for t=1:n
Bnew(2:r-1,2:c-1) =K(1)*Bold(1:r-2,1:c-2)+K(2)*Bold(1:r-2,2:c-1)+K(3)*Bold(1:r-2,3:c);
Bnew(2:r-1,2:c-1)+=K(4)*Bold(2:r-1,1:c-2)+K(5)*Bold(2:r-1,2:c-1)+K(6)*Bold(2:r-1,3:c);
Bnew(2:r-1,2:c-1)+=K(7)*Bold(3:r,1:c-2)+K(8)*Bold(3:r,2:c-1)+K(9)*Bold(3:r,3:c);
% on se ramene a [0;255], inutile de diviser avant
m=min(Bnew(:));
M=max(Bnew(:));
Bnew=255/(M-m).*(Bnew-m);
Bold=Bnew;
end
end

Puis on fait appelle à cette fonction en changeant le masque et le nombre d’itérations :

1. clear all B = myflou(A

,[1/9,1/9,1/9,1/9,1/9,1/9,1/9,1/9,1/9],50);
A = double(imread('lena512.bmp'));
% on se ramene a [0;255] subplot(1,3,1)
m = min(A(:)); imshow(uint8(A));
M = max(A(:)); title ( "Original" );
A = 255/(M-m).*(A-m); %imwrite(uint8(A),'exo5orig.jpg','jpg');
colormap(gray(256));
subplot(1,3,2)

© 2022-2023 G. Faccanoni 59
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

imshow(uint8(B)); subplot(1,3,3)
title ( "Moyenne 9" ); imshow(uint8(abs(B-A)));
%imwrite(uint8(B),'exo5flout.jpg','jpg'); title ( "|moy(A)-A|" );
%imwrite(uint8(abs(B-A)),'exo5err.jpg','jpg');

2. clear all subplot(1,3,1)

imshow(uint8(A));
A = double(imread('lena512.bmp')); title ( "Original" );
% on se ramene a [0;255] %imwrite(uint8(A),'exo52orig.jpg','jpg');
m = min(A(:));
M = max(A(:)); subplot(1,3,2)
A = 255/(M-m).*(A-m); imshow(uint8(B));
colormap(gray(256)); title ( "Filtre rehausseur de contraste" );
%imwrite(uint8(B),'exo52flout.jpg','jpg');
%B = myflou(A
,[-1/9,-1/9,-1/9,-1/9,1,-1/9,-1/9,-1/9,-1/9],1)subplot(1,3,3)
; imshow(uint8(abs(B-A)));
B = myflou(A,[0,-1/2,0,-1/2,3,-1/2,0,-1/2,0],1); title ( "|reh(A)-A|" );
%imwrite(uint8(abs(B-A)),'exo52err.jpg','jpg');

3. clear all title ( "Original" );

%imwrite(uint8(A),'exo53orig.jpg','jpg');
A = double(imread('lena512.bmp'));
% on se ramene a [0;255] subplot(1,3,2)
m = min(A(:)); imshow(uint8(B));
M = max(A(:)); title ( "Detection des contours" );
A = 255/(M-m).*(A-m); %imwrite(uint8(B),'exo53flout.jpg','jpg');
colormap(gray(256));
subplot(1,3,3)
B = myflou(A,[-1,-1,-1,-1,8,-1,-1,-1,-1],1); imshow(uint8(abs(B-A)));
title ( "|dc(A)-A|" );
subplot(1,3,1) %imwrite(uint8(abs(B-A)),'exo53err.jpg','jpg');
imshow(uint8(A));

4. clear all
subplot(1,2,1)
A = double(imread('lena512.bmp')); imshow(uint8(A));
% on se ramene a [0;255] title ( "Original" );
m = min(A(:)); %imwrite(uint8(A),'5rehaussement.jpg','jpg');
M = max(A(:));
A = 255/(M-m).*(A-m); subplot(1,2,2)
colormap(gray(256)); N = A-0.1*B;
imshow(uint8(N));
B = myflou(A title ( "Rheaussement" );
,[1/8,1/8,1/8,1/8,-1,1/8,1/8,1/8,1/8],1); %imwrite(uint8(N),'5rehaussement.jpg','jpg');

5. clear all ,[1/8,1/4,1/8,0,0,0,-1/8,-1/4,-1/8],1);

G = sqrt(B1.^2+B2.^2);
A = double(imread('lena512.bmp'));
% on se ramene a [0;255] subplot(1,2,1)
m = min(A(:)); imshow(uint8(A));
M = max(A(:)); title ( "Original" );
A = 255/(M-m).*(A-m);
colormap(gray(256)); subplot(1,2,2)
imshow(uint8(G));
B1 = myflou(A title ( "Contours" );
,[1/8,0,-1/8,1/4,0,-1/4,1/8,0,-1/8],1); %imwrite(uint8(G),'6cont.jpg','jpg');
B2 = myflou(A

60 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.6 Images à couleurs

2.6 Images à couleurs

Une image couleur est en réalité composée de trois images, afin de représenter le rouge, le vert, et le bleu. Chacune de
ces trois images s’appelle un canal. Cette représentation en rouge, vert et bleu mime le fonctionnement du système visuel
humain.
La figure suivante montre la décomposition d’une image couleur en ses trois canaux constitutifs.

(a) Originale (b) Canal rouge (c) Canal Vert (d) Canal Bleu

F IGURE 2.13 – Canaux RGB

Chaque pixel de l’image couleur contient ainsi trois nombres (r, g , b), chacun étant un nombre entier entre 0 et 255. Un tel
nombre est représentable sur 8 bits et s’appelle un octet. Il y a donc 2563 = 224 = 16777216 couleurs possibles.
Si le pixel est égal à (r, g , b) = (255, 0, 0), il ne contient que de l’information rouge, et est affiché comme du rouge. De façon
similaire, les pixels valant (0, 255, 0) et (0, 0, 255) sont respectivement affichés vert et bleu.

On peut afficher à l’écran une image couleur à partir de ses

trois canaux (r, g , b) en utilisant les règles de la synthèse ad-
ditive des couleurs : le triplet (0, 0, 0) correspond à un pixel
noir alors qu’un pixel blanc est donné par (255, 255, 255).
La figure suivante montre les règles de composition cette
synthèse additive des couleurs. Un pixel avec les valeurs
(r, g , b) = (255, 0, 255) est un mélange de rouge et de vert, il est
ainsi affiché comme jaune. F IGURE 2.14 – Synthèse additive des couleurs

Une autre représentation courante pour les images couleurs utilise comme couleurs de base le cyan, le magenta et le jaune.
On calcule les trois nombres (c, m, j ) correspondant à chacun de ces trois canaux à partir des canaux rouge, vert et bleu
(r, v, b) comme suit 
 c = 255 − r,

m = 255 − v,

j = 255 − b.


Par exemple, un pixel de bleu pur (r, v, b) = (0, 0, 255) va devenir (c, m, j ) = (255, 255, 0). La figure suivante montre les trois
canaux (c, m, j ) d’une image couleur.

(a) Originale (b) Canal cyan (c) Canal magenta (d) Canal jaune

F IGURE 2.15 – Canaux CMY

Afin d’afficher une image couleur à l’écran à partir des trois

canaux (c, m, j ), on doit utiliser la synthèse soustractive des
couleurs. La figure suivante montre les règles de composi-
tion cette synthèse soustractive. Notons que ces règles sont
celles que l’on utilise en peinture, lorsque l’on mélange des
pigments colorés. Le cyan, le magenta et le jaune sont appelés
couleurs primaires.
F IGURE 2.16 – Synthèse soustractive des couleurs

© 2022-2023 G. Faccanoni 61
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

On peut donc stocker sur un disque dur une image couleur en stockant les trois canaux, correspondant aux valeurs (r, g , b)
ou (c, m, j ).
On peut modifier les images couleur tout comme les images en niveaux de gris. La façon la plus simple de procéder consiste
à appliquer la modification à chacun des canaux.

Exercice 2.20 (Devine le résultat)

Devine le résultat :
clear all
mat = zeros(2,2,3); % images en couleurs
mat(1,1,1) = 255; % mat(i,j,1) = > R(i,j)
mat(2,1,2) = 255; % mat(i,j,2) = > G(i,j)
mat(1,2,3) = 255; % mat(i,j,3) = > B(i,j)
%mat(2,2,:) = 255;
%mat(2,2,[1,2]) = 255;
%mat(2,2,[1,3]) = 255;
%mat(2,2,[2,3]) = 255;
imshow(uint8(mat));

Correction
Il s’agit des trois matrices R, G, B suivantes (correspondantes aux trois canaux RGB) :
µ ¶ µ ¶ µ ¶
255 0 0 0 0 255
R= G= B=
0 0 255 0 0 0

ainsi

255r + 0g + 0b 0r + 0g + 255b red blue

=
0r + 255g + 0b 0r + 0g + 0b green black

Exercice 2.21 (Affichage image couleur)

Soit les trois matrices
     
0 64 128 64 128 192 128 192 255
R =  64 128 192 G = 128 192 255 B= 0 64 128
128 192 255 0 64 128 64 128 192

Après avoir défini l’hypermatrice A par l’instruction A=cat(3, R,G,B ); afficher l’image associée.

Correction

clear all

R = [0 64 128; 64 128 192; 128 192 255];

G = [64 128 192; 128 192 255; 0 64 128];
B = [128 192 255; 0 64 128; 64 128 192];
A = cat(3, R,G,B );
imshow(uint8(A));
%imwrite(uint8(A),'exo1Color.jpg','jpg');

Exercice 2.22
On considère l’image 2.17a. Obtenir les autres images par permutations des canaux.

62 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.6 Images à couleurs

(a) RGB (b) BRG (c) GRB (d) RBG (e) BGR (f) GBR

F IGURE 2.17 – Permutations

Correction N = cat(3, G,R,B );

clear all imshow(uint8(N));
title ( "GRB" );
A = double(imread('hibiscus.png')); %imwrite(uint8(N),'exo2ColorGRB.jpg','jpg');
% size(A)
R = A(:,:,1); subplot(2,3,4)
G = A(:,:,2); N = cat(3, R,B,G );
B = A(:,:,3); imshow(uint8(N));
title ( "RBG" );
subplot(2,3,1) %imwrite(uint8(N),'exo2ColorRBG.jpg','jpg');
imshow(uint8(A));
title ( "RGB" ); subplot(2,3,5)
%imwrite(uint8(A),'exo2ColorRGB.jpg','jpg'); N = cat(3, B,G,R );
imshow(uint8(N));
subplot(2,3,2) title ( "BGR" );
N = cat(3, B,R,G ); %imwrite(uint8(N),'exo2ColorBGR.jpg','jpg');
% ou, en modifiant A
% A(:,:,[1,2,3]) = A(:,:,[3,1,2]); subplot(2,3,6)
imshow(uint8(N)); N = cat(3, G,B,R );
title ( "BRG" ); imshow(uint8(N));
%imwrite(uint8(N),'exo2ColorBRG.jpg','jpg'); title ( "GBR" );
%imwrite(uint8(N),'exo2ColorGBR.jpg','jpg');
subplot(2,3,3)

Exercice 2.23 (Négatif )

On considère l’image 2.18a. Obtenir l’image 2.18b en considérant le négatif pour tous les canaux.

(a) Original (b) Canaux néga-

tifs

F IGURE 2.18 – Négatif

Correction %imwrite(uint8(A),'exo5ColorRGB.jpg','jpg');
clear all
subplot(1,2,2)
A = double(imread('hibiscus.png')); N = 255-A;
size(A) imshow(uint8(N));
title ( "Negatif" );
subplot(1,2,1) %imwrite(uint8(N),'exo5ColorRGBNeg.jpg','jpg');
imshow(uint8(A));
title ( "RGB" );

© 2022-2023 G. Faccanoni 63
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

Exercice 2.24 (Modification d’un seul canal)

On considère l’image 2.19a. Obtenir l’image 2.19b en considérant le négatif pour le canal rouge.

(a) Original (b) Canal rouge

négatif

F IGURE 2.19 – Modification du canal de rouge

Correction %imwrite(uint8(A),'exo3ColorRGB.jpg','jpg');
clear all
subplot(1,2,2)
A = double(imread('hibiscus.png')); M = (255-A(:,:,1));
size(A) N = cat(3, M,A(:,:,2),A(:,:,3));
imshow(uint8(N));
subplot(1,2,1) title ( "Modification canal rouge" );
imshow(uint8(A)); %imwrite(uint8(N),'exo3Color0GB.jpg','jpg');
title ( "RGB" );

On peut calculer une image en niveaux de gris à partir d’une

image couleur en moyennant les trois canaux :

r i j + g i j + bi j
bi j =
3
qui s’appelle la luminance de la couleur. La figure ci-contre
(a) Originale (b) Luminance
montre l’image de luminance associée à une image couleur.
F IGURE 2.20 – Luminance

Exercice 2.25 (RGB → niveau de gris)

On considère l’image RGB 2.21a.
① Obtenir l’image en niveaux de gris 2.21b en considérant, pour chaque pixel, la moyenne arithmétique des trois
canaux (luminance).
② L’œil est plus sensible à certaines couleurs qu’à d’autres. Le vert (pur), par exemple, paraît plus clair que le bleu
(pur). Pour tenir compte de cette sensibilité dans la transformation d’une image couleur en une image en niveaux
de gris, on ne prend généralement pas la moyenne arithmétique des intensités de couleurs fondamentales, mais
une moyenne pondérée. La formule standard donnant le niveau de gris de chaque pixel en fonction des trois
composantes est
gray = E (0.21 · r + 0.72 · g + 0.07 · b)
où E (·) désigne la partie entière. Obtenir l’image en niveaux de gris 2.21c selon cette moyenne.
③ Une autre moyenne consiste à prendre, pour chaque pixel, la moyenne du canal le plus présent et celui le moins
présent : © ª © ª¶
max r, g , b + min r, g , b
µ
gray = E .
2
Obtenir l’image en niveaux de gris 2.21d selon cette moyenne.

64 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.6 Images à couleurs

(a) Original (b) Moyenne arithmé- (c) Moyenne pondé- (d) Moyenne max-min
tique rée

F IGURE 2.21 – RGB → niveau de gris

Correction imshow(uint8(N));
clear all title ( "Aritmetique" );
%imwrite(uint8(N),'exo6Color1.jpg','jpg');
A = double(imread('hibiscus.png'));
[row,col,couches] = size(A) subplot(1,4,3)
N = floor(0.21*R+0.72*G+0.07*B);
R = A(:,:,1); imshow(uint8(N));
G = A(:,:,2); title ( "Ponderee" );
B = A(:,:,3); %imwrite(uint8(N),'exo6Color2.jpg','jpg');

subplot(1,4,1) subplot(1,4,4)
imshow(uint8(A)); N = ( max(max(R,G),B) + min(min(R,G),B) ) /2 ;
title ( "RGB" ); imshow(uint8(N));
%imwrite(uint8(A),'exo6Color.jpg','jpg'); title ( "Max-Min" );
%imwrite(uint8(N),'exo6Color3.jpg','jpg');
subplot(1,4,2)
N = floor(R+G+B)/3;

Exercice 2.26 (Fusion de deux images)

Considérons les deux images 2.22a et 2.22e (qui ont la même dimension). Superposer les deux images. Créer un film qui
montre le passage de l’une à l’autre.

(a) Première (b) Intermédiaire (c) Intermédiaire (d) Intermédiaire (e) Deuxième
image image

F IGURE 2.22 – Morphing

Correction morphing = (k/N)*im1+(1-k/N)*im2; % NE PAS

clear all; clc; APPELER LE FICHIER COMME LA MATRICE!
imshow(morphing);
[fig1,map1] = imread("s1.png"); pause(0.5);
im1 = ind2rgb(fig1,map1); end
[fig2,map2] = imread("s2.png"); %imwrite(0.9*im1+0.1*im2,'s12a.png','png')
im2 = ind2rgb(fig2,map2); %imwrite(0.5*im1+0.5*im2,'s12.png','png')
N = 16; %imwrite(0.1*im1+0.9*im2,'s12c.png','png')
for k = 1:N
%subplot(4,4,k)

© 2022-2023 G. Faccanoni 65
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

2.7 Décomposition en valeurs singulières et compression JPG

Tout comme pour la réduction du nombre de pixels, la réduction du nombre de niveaux de gris influe beaucoup sur la qualité
de l’image. Afin de réduire au maximum la taille d’une image sans modifier sa qualité, on utilise des méthodes plus complexes
de compression d’image. Une des première méthode s’appelle JPEG, basée sur la décomposition en valeurs singulières. La
méthode la plus efficace s’appelle JPEG-2000. Elle utilise la théorie des ondelettes. Pour en savoir plus à ce sujet, vous pouvez
consulter cet article d’Erwan Le Pennec.
Soit A ∈ Rn×p une matrice rectangulaire. Un théorème démontré officiellement en 1936 par C. E CKART et G. Y OUNG affirme
que toute matrice rectangulaire A se décompose sous la forme

A = USVT

avec U ∈ Rn×n et V ∈ Rp×p des matrices orthogonales T T

ª (i.e. U = U et V = V ) et S ∈ R
−1 −1 n×p
une matrice diagonale qui
contient les r valeurs singulières de A, r = min n, p , σ1 ≥ σ2 ≥ σr ≥ 0. Ce qui est remarquable, c’est que n’importe quelle
©

matrice admet une telle décomposition alors que la décomposition en valeurs propres (la diagonalisation d’une matrice)
n’est pas toujours possible.
Notons ui et vi les vecteurs colonne des matrices U et V. La décomposition s’écrit alors

σ1 vT1
  

..  . 
  .. 


 .  
  vT 
¢ σr
A = USVT = u1
¡  r 
... ur ur +1 ... un 

 T 
| {z } 0  vr +1 
n×n ..   .. 
  

 .  . 
0 vTp
| {z } | {z }
n×p p×p
 T 
σ1 v1

r
..   ..  X
σi ui × vTi
¡ ¢
= u1 ... ur  .  .  =
i =1
σr vTr
| {z } | {z }
n×r r ×r
| {z } | {z }
r ×r r ×p

Pour calculer ces trois matrices on remarque que A = USVT = USV−1 et AT = VSUT = VSU−1 ainsi, pour i = 1, . . . , r , en
multipliant par A à gauche AT ui = σi vi et en multipliant par A à droite Avi = σi ui on obtient

AAT ui = σi Avi = σ2i ui , pour i = 1, . . . , r

AT Avi = σi AT ui = σ2i vi , pour i = 1, . . . , r

ainsi les σ2i sont les valeurs propres de la matrice AAT et les ui les vecteurs propres associés mais aussi les σ2i sont les valeurs
propres de la matrice AT A et les vi les vecteurs propres associés (attention, étant des valeurs propres, ils ne sont pas définis
de façon unique).
On peut exploiter cette décomposition pour faire des économies de mémoire.
⋆ Pour stocker la matrice A nous avons besoin de n × p valeurs.
⋆ Pour stocker la décomposition SVD nous avons besoin de n × r + r + r × p = (n + p + 1)r > (n + p + 1)r valeurs donc à
priori on ne fait pas d’économies de stockage. Cependant, s’il existe s < r tel que σs = σs+1 = · · · = σr = 0, alors nous
n’avons plus besoin que de n × s + s + s × p = (n + p + 1)s valeurs. Si s < np/(n + p + 1) on fait des économies de stockage.
Idée de la compression : si nous approchons A en ne gardant que les premiers s termes de la somme (sachant que les derniers
termes sont multipliés par des σi plus petits, voire nuls)
s
A σi ui × vTi ,
X
e= où s < r
i =1 | {z }
∈Rn×p

⋆ pour stocker la matrice A

e nous avons toujours besoin de n × p valeurs,
⋆ pour stocker la décomposition SVD nous avons besoin de n × s + s + s × p = (n + p + 1)s valeurs. Si s < np/(n + p + 1) on
fait des économies de stockage.

66 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.7 Décomposition en valeurs singulières et compression JPG

Exercice 2.27 (Valeurs singulières)

Calculer analytiquement et vérifier numériquement sa décomposition SVD de la matrice
µ ¶
1 2 0
A=
2 1 0

Correction
A ∈ Rn×p avec n = 2 et p = 3 donc r = 2.
Pour calculer la décomposition SVD nous allons calculer les valeurs et vecteurs propres des matrices AAT et AT A.

Valeurs propres Vecteurs propres unitaires

µ ¶ µ ¶
T 5 4 1 1 1
AA = λ1 = 9 > λ2 = 1 U= p
4 5 2 1 −1
1 1 0
   
5 4 0
T 1 
A A = 4 5 0  λ1 = 2 > λ2 = 1 > λ3 = 0 V = p 1 −1 p0 
0 0 0 2 0 0 2

Donc

σ1 vT1
  

..  . 
  .. 


 .  
  vT 
¢ σr
A = USVT = u1
¡  r 
... ur ur +1 ... un 

 T 
| {z } 0  vr +1 
∈Rn×n ..   .. 
  

 .  . 
0 vTp
| {z } | {z }
∈Rn×p ∈Rp×p
σ1 vT1
  
r
..   ..  X
σi ui × vTi
¡ ¢
= u1 ... ur  .  .  =
i =1
σr vTr
| {z } | {z }
∈Rn×r ∈Rr ×r
| {z } | {z }
∈Rr ×r ∈Rr ×p

devient

1 1 0
 
µ ¶µ ¶
1 1 1 3 0 0 1 
A= p p 1 −1
p
0 
2 1 −1 0 1 0 2 0 0 2
µ ¶µ ¶µ ¶
1
r =2 1 1 3 0 1 1 0
=
2 1 −1 0 1 1 −1 0
µ ¶ µ ¶
3 1 1 0 1 1 −1 0
= +
2 1 1 0 2 −1 1 0

Notons que la décomposition n’est pas unique, par exemple avec Octave on trouve

Vecteurs propres unitaires :

1 0
 
µ ¶ −1
1 1 −1 1 
U= p V= p 1 1
p
0 
2 1 1 2 0 0 2

ce qui donne le même résultat (heureusement !)

A=[1 2 0; 2 1 0]
[n,p]=size(A)
r=min(n,p)

AAT=A*A'
[VecAAT,ValAAT]=eig(AAT) % unsorted list of all eigenvalues
% To produce a sorted vector with the eigenvalues, and re-order the eigenvectors accordingly:

© 2022-2023 G. Faccanoni 67
Chapitre 2 Traitement mathématique des images numériques Mis à jour le Lundi 28 août 2023

[ee,perm] = sort(diag(abs(ValAAT)),"descend");
ValAAT=diag(ee)
VecAAT=VecAAT(:,perm)

ATA=A'*A
[VecATA,ValATA]=eig(ATA)
[ee,perm] = sort(diag(abs(ValATA)),"descend");
ValATA=diag(ee)
VecATA=VecATA(:,perm)

myS=diag(sqrt(diag(ValATA)),n,p)
myU=VecAAT
myV=VecATA

[UU,SS,VV]=svd(A)

dS=diag(SS)

AA=zeros(5,4);
for i=1:numel(dS)
temp=dS(i)*UU(:,i)*VV(i,:)
AA+=temp
end

Exercice 2.28 (Compression par SVD)

1. Tester la compression avec s = 10 et s = 100 pour obtenir les images 2.23 ainsi que la carte des erreurs.
n ¯ o
np
2. Calculer max s ∈ [0; r ] ¯ s < n+p+1 qui est la limite en dessous de laquelle le stockage de la décomposition SVD
¯

permet des économies de mémoire par rapport au stockage de l’image initiale.

3. Supposons que la précision d’Octave soit de l’ordre de 3 chiffres significatifs. Alors les valeurs singulières signi-
ficatives doivent avoir un rapport de moins de 10−3 avec la valeur maximale σ1 , les autres étant considérées
comme «erronées». Calculer le nombre de valeurs singulières «significatives», i.e. la plus grande valeur de s telle
σ
que σ1i < 10−3 pour i = s + 1, . . . , r .

Correction
clear all n = 2;
for s = [vsSignif,floor(economie),100,10]
A = double(imread('lena512.bmp')); subplot(5,2,2*n-1)
[row,col] = size(A) X = U(:,1:s)*S(1:s,1:s)*(V(:,1:s))';
imshow(uint8(X));
colormap(gray(256)); title ( strcat("s = ",num2str(s)) );
%imwrite(uint8(X),strcat(['exo6-' num2str(s) '.
subplot(5,2,1) jpg']),'jpg');
imshow(uint8(A)); subplot(5,2,2*n)
s = min(row,col); erreur = abs(X-A);
title ( strcat("s = ",num2str(s)) ); somerr = sum(erreur(:));
imwrite(uint8(A),strcat(['exo6-' num2str(s) '.jpg' m = min(erreur(:));
]),'jpg'); M = max(erreur(:));
erreur = 255-255/(M-m)*(erreur-m);
[U,S,V] = svd(A); imshow(uint8(erreur));
%imwrite(uint8(erreur),strcat(['exo6-' num2str(
% on fait des economies de stockage si "s" est < a s) 'E.jpg']),'jpg');
"economie" : n+ = 1;
economie = row*col/(row+col+1) end
vsSignif = sum(sum( (S./S(1,1))>1.e-3 ))
n ¯ o
np
On a max s ∈ [0; r ] ¯ s < n+p+1 = 255 : on fait des économies de stockage tant qu’on garde au plus les premières 255 valeurs
¯

singulières.
La photo de Lena, de taille 512 × 512, possède 275 valeurs singulières «significatives».

68 © 2022-2023 G. Faccanoni
Mis à jour le Lundi 28 août 2023 2.7 Décomposition en valeurs singulières et compression JPG

(a) Original s = n = p = 512

¯σ
¯ ¯
np
n o n o
(b) min s ∈ [0; r ] ¯ σ s < 10−5 = (c) max s ∈ [0; r ] ¯ s < n+p+1 =
¯ (d) s = 100 (e) s = 10
1
275 255

¯σ
¯ ¯
np
n o n o
(f) min s ∈ [0; r ] ¯ σ s < 10−5 = (g) max s ∈ [0; r ] ¯ s < n+p+1 =
¯ (h) s = 100 (i) s = 10
1
275 255

F IGURE 2.23 – SVD - exercice 2.28

© 2022-2023 G. Faccanoni 69