Traitement d'images avec R

application du calcul matriciel

1 Introduction

L’objectif de cette activité est d’illustrer le rôle du calcul matriciel (applications linéaires, transformations) et des outils géométriques (produit scalaire, norme) en traitement numérique des images. Chaque exercice consiste à créer un filtre modifiant le contenu d’un fichier image.

Vous pouvez télécharger les deux fichiers PNG ci-dessous, utilisés dans les exemples, pour tester l’ensemble de vos filtres :

2 Le package PNG

Ce package R permet de charger des fichiers image au format PNG (Portable Network Graphics). La fonction readPNG initialise une variable de type array de taille hauteur × largeur × profondeur : hauteur et largeur correspondent aux dimensions, en pixels, de l’image chargée. La profondeur, quant à elle, correspond au nombre de composantes du pixel :

Chaque composante est un nombre réel compris dans l’intervalle [0,1].

Pour l’utiliser, vous devez cocher le package PNG dans la liste des packages de RStudio ou l’importer à l’aide de la fonction library.

# Importation du package PNG
library("png")
# Exemple de chargement d'un fichier PNG
# localisé dans le sous-répertoire img
img1<-readPNG("img/solo-256px.png")

2.1 Dimensions de l’image

La fonction ncol permet d’obtenir la largeur de l’image exprimée en pixels.

width<-ncol(img1)
print(paste("width: ",width))
## [1] "width:  256"

La fonction nrow permet d’obtenir la hauteur de l’image exprimée en pixels.

height<-nrow(img1)
print(paste("height: ",height))
## [1] "height:  256"

2.2 Accès à un pixel de l’image

L’image étant stockée en mémoire sous la forme d’un objet array, l’accès à la composante numéro \(c\) du pixel de coordonnées \((x,y)\) se fait via l’opérateur [y,x,c]. Par exemple, le code ci-dessous affiche les composantes rouge, vert, bleu du pixel de coordonnées \(x=42\) et \(y=57\).

print(img1[57,42,])
## [1] 1.0000000 0.9058824 0.4235294

Le pixel du coin supérieur gauche de l’image a pour coordonnées \((1,1)\), l’axe \(x\) est horizontal dirigé vers la droite et l’axe \(y\) est vertical dirigé vers le bas. Ainsi le pixel du coin inférieur droit de l’image a pour coordonnées \((w,h)\).

2.3 Affichage de l’image

Pour cela, vous pouvez utiliser la fonction display(...) développée par Frédéric Blanchard qui permet d’afficher plusieurs images dans un seul graphique. Cette fonction est disponible dans le script color_utils.R.

# Affichage de l'image img1 
display(img1)

3 Modification de couleurs

3.1 Niveaux de gris

Pour convertir une image couleur RGB en une image en niveaux de gris, il faut calculer la luminance de chaque pixel. Cette luminance est une combinaison linéaire des composantes rouge, vert et bleu :

\[ \left( \begin{array}{c} r' \\ g' \\ b' \end{array} \right) = \begin{pmatrix} 0.2126 & 0.7152 & 0.0722 \\ 0.2126 & 0.7152 & 0.0722 \\ 0.2126 & 0.7152 & 0.0722 \end{pmatrix} \left( \begin{array}{c} r \\ g \\ b \end{array} \right) \]

Exercice 1
Écrire une fonction R greyscale(src) qui retourne une image en niveaux de gris à partir de l’image RGB src passée en paramètre.

display(img1,greyscale(img1))

3.2 Sépia

Le ton sépia rappelle l’aspect visuel des photographies du début du XXe siècle. Cet effet peut être obtenu grâce à l’application linéaire dont la matrice figure ci-dessous :

\[ \left( \begin{array}{c} r' \\ g' \\ b' \end{array} \right) = \begin{pmatrix} 0.393 & 0.769 & 0.189 \\ 0.349 & 0.686 & 0.168 \\ 0.272 & 0.534 & 0.131 \end{pmatrix} \left( \begin{array}{c} r \\ g \\ b \end{array} \right) \]

Dans la mesure où la somme des coefficients d’une même ligne est supérieure à 1, le résultat de la combinaison linéaire doit être limité à 1 si celui-ci est hors de l’intervalle \([0,1]\). Pour cela, vous pouvez utiliser la fonction clamp(pixel) disponible dans le script color_utils.R.

#
# Limit pixel values in the [0,1] interval
#
# Parameter:
# - pixel: array containing real numbers
#
# Returns:
#  array containing real numbers between 0 and 1
#
clamp<-function(pixel){
  for(i in 1:length(pixel)){
    if(pixel[i]<0){
      pixel[i]<-0
    }
    else if(pixel[i]>1){
      pixel[i]<-1
    }
  }
  return(pixel)
}

Exercice 2
Écrire une fonction R sepia(src) qui retourne une image au ton sépia à partir de l’image RGB src passée en paramètre.

display(img1,sepia(img1))

3.3 Mélange d’images

Cette dernière utilisation des combinaisons linéaires consiste à mélanger deux images de même taille en calculant une moyenne pondérée des pixels des deux images source. La formule ci-dessous illustre le calcul à réaliser :

\[ \left( \begin{array}{c} r' \\ g' \\ b' \end{array} \right) = \alpha \left( \begin{array}{c} r_1 \\ g_1 \\ b_1 \end{array} \right) + (1-\alpha) \left( \begin{array}{c} r_2 \\ g_2 \\ b_2 \end{array} \right) \quad \text{où} \quad \alpha \in [0,1] \]

Exercice 3
Écrire une fonction R mix(src1,src2,factor) qui retourne une image issue du mélange des images src1 et src2 passées en paramètres avec un poids égal à factor.

img2<-readPNG("img/logo-starwars-256px.png")
img3<-mix(img1,img2,0.5)
display(img1,img3,img2)

Une autre manière de mélanger des images consiste à multiplier les composantes des pixels d’une image avec celles d’un masque (i.e. une image de luminance). La formule ci-dessous illustre le calcul à réaliser :

\[ \left( \begin{array}{c} r' \\ g' \\ b' \end{array} \right) = \left( \begin{array}{c} r \cdot r_m \\ g \cdot g_m \\ b \cdot g_m \end{array} \right) \]

Exercice 3 bis
Écrire une fonction R multiply(src,mask) qui retourne une image issue du produit pixel par pixel d’une image src avec un masque mask.

img3<-multiply(img1,img2)
display(img1,img3,img2)

3.4 Convolution

De façon très concrète, la convolution consiste à calculer la couleur d’un pixel de l’image destination comme la moyenne pondérée du pixel de l’image source et de ses voisins. En général, le voisinage est carré, centré sur le pixel courant. Les poids utilisés pour calculer la moyenne sont stockés dans une matrice appelée masque de convolution.

Exercice 4
Écrire une fonction R convolution(src,mask) qui convolue une image src passée en paramètre en utilisant un masque de convolution mask représenté par une matrice 3x3. Pour des considérations de simplicité, on ne traitera pas les pixels se trouvant sur la bordure de l’image destination

\[ \text{Flou moyen} \qquad \textbf{M}_1 = \begin{pmatrix} 1/9 & 1/9 & 1/9 \\ 1/9 & 1/9 & 1/9 \\ 1/9 & 1/9 & 1/9 \end{pmatrix} \]

mask<-matrix(c(1/9,1/9,1/9,1/9,1/9,1/9,1/9,1/9,1/9),nrow=3,byrow=TRUE)
display(img1,convolution(img1,mask))

\[ \text{Réhaussement de contraste} \qquad \textbf{M}_2 = \begin{pmatrix} -1 & -1 & -1 \\ -1 & 9 & -1 \\ -1 & -1 & -1 \end{pmatrix} \]

mask<-matrix(c(-1,-1,-1,-1,9,-1,-1,-1,-1),nrow=3,byrow=TRUE)
display(img1,convolution(img1,mask))

\[ \text{Détection de contours (Sobel)} \qquad \textbf{M}_3 = \begin{pmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{pmatrix} \]

mask<-matrix(c(-1,0,1,-2,0,2,-1,0,1),nrow=3,byrow=TRUE)
display(img1,convolution(img1,mask))

Pour certains filtres de convolution, il peut être intéressant de travailler avec des masques de plus grande taille. C’est le cas du filtre flou moyen par exemple. Plus le voisinage est grand, plus l’image est floue.

Exercice 5
Écrire une fonction R blur(src,halfw) qui floute une image src passée en paramètre en utilisant un voisinage carré de demi-largeur halfw.

display(img1,blur(img1,3))

Le filtre flou peut devenir lent avec des tailles de voisinage importantes. Ce filtre étant séparable, il peut se décomposer en deux filtres consécutifs : un flou horizontal suivi d’un flou vertical.

Exercice 6
Écrire deux fonctions R hblur(src,halfw) et vblur(src,halfw) qui floutent une image src passée en paramètre en parcourant un intervalle horizontal (resp. vertical) de demi-largeur halfw.

display(img1,hblur(img1,7),vblur(img1,7))

Le graphique ci-dessous représente l’évolution des temps de calcul du filtre flou blur (courbe rouge) et de la combinaison du filtre flou horizontal hblur avec le filtre flou vertical vblur (courbe verte).

On définit la fonction smartblur(src,halfw) qui applique successivement un flou vertical et un flou horizontal.

smartblur<-function(src,halfw){
  return(hblur(vblur(src,halfw),halfw))
}

display(blur(img1,7),smartblur(img1,7))

3.4.1 Application au débruitage d’images

Le filtre flou peut être utilisé pour atténuer le bruit au sein d’une image en remplaçant la valeur d’un pixel par la moyenne de ses voisins. Néanmoins cette technique s’avère inefficace avec un bruit ponctuel sel et poivre.

saltnpepper<-function(src,density){
  w<-ncol(src)
  h<-nrow(src)
  n<-floor(w*h*density/2)

  for(i in 1:n){
    x<-runif(1,1,w)
    y<-runif(1,1,h)
    src[x,y,]<-c(0,0,0)

    x<-runif(1,1,w)
    y<-runif(1,1,h)
    src[x,y,]<-c(1,1,1)
  }

  return(src)
}

img1n<-saltnpepper(img1,0.1)
display(img1n,smartblur(img1n,1),smartblur(img1n,3))

Le filtre médian quant à lui consiste à rassembler dans un tableau les valeurs des pixels voisins du pixel courant. La valeur du pixel résultat correspond alors à la médiane de ce tableau.

Pour des images RGB, le filtre médian est appliqué indépendamment sur chacune des 3 composantes à l’aide de la fonction R median().

Exercice 7
Écrire une fonction R medianFilter(src,halfw) qui applique un filtre médian avec un voisinage carré centré de demi-largeur halfw à une image src passée en paramètre.

img1n<-saltnpepper(img1,0.1)
display(img1n,medianFilter(img1n,1))

3.5 Modèles de couleur : RGB, HSL

3.5.1 Modèle RGB (Red Green Blue)

Ce modèle repose sur le principe de la synthèse additive qui permet de générer une grande variété de couleurs à partir des trois couleurs rouge, vert et bleu. En utilisant 256 nuances pour chaque couleur, on peut ainsi produire jusqu’à 2563 = 16 777 216 couleurs différentes.

Représentation graphique de l’espace du modèle RGB

Représentation graphique de l’espace du modèle RGB

3.5.2 Modèle HSL (Hue Saturation Lightness)

HSL est un modèle de couleur alternatif à RGB. La teinte (hue) est représentée par un angle exprimé en degrés couvrant tout le spectre colorimétrique. La saturation, exprimée en pourcentage, permet de jouer sur le contraste, enfin la luminance (lightness) permet de définir une couleur plus ou moins claire.

Représentation graphique de l’espace du modèle HSL

Représentation graphique de l’espace du modèle HSL

3.5.3 Fonctions de conversion

Afin de faciliter le passage d’un modèle de couleur à l’autre, vous pourrez utiliser les fonctions utilitaires rgb2hsl(pixel) et hsl2rgb(pixel) dont le code est disponible dans le script color_utils.R.

#
# RGB to HSL conversion function
# http://www.rapidtables.com/convert/color/rgb-to-hsl.htm
#
# Parameter:
# - pixel: array containing red, green and blue values
#
# Returns:
#  array containing hue, saturation and lightness values
#
rgb2hsl<-function(pixel){
  r<-pixel[1]
  g<-pixel[2]
  b<-pixel[3]

  cmax<-max(r,g,b)
  cmin<-min(r,g,b)
  delta<-cmax-cmin

  # Lightness calculation
  l<-(cmax+cmin)/2

  # Saturation calculation
  s<-0
  if(delta!=0){
    s<-delta/(1-abs(2*l-1))
  }

  # Hue calculation
  h<-0
  if(delta!=0){
    if(cmax==r){
      h<-(((g-b)/delta)%%6)*60
    }
    if(cmax==g){
      h<-(2+(b-r)/delta)*60
    }
    if(cmax==b){
      h<-(4+(r-g)/delta)*60
    }
  }
  return(c(h,s,l))
}

#
# HSL to RGB conversion function
# http://www.rapidtables.com/convert/color/hsl-to-rgb.htm
#
# Parameter:
# - pixel: array containing hue, saturation and lightness values
#
# Returns:
#  array containing red, green and blue values
#
hsl2rgb<-function(pixel){
  h<-pixel[1]
  s<-pixel[2]
  l<-pixel[3]

  c<-(1-abs(2*l-1))*s
  x<-c*(1-abs((h/60)%%2-1))

  r<-0
  g<-0
  b<-0
  if((h>=0)&(h<60)){
    r<-c
    g<-x
  }
  if((h>=60)&(h<120)){
    r<-x
    g<-c
  }
  if((h>=120)&(h<180)){
    g<-c
    b<-x
  }
  if((h>=180)&(h<240)){
    g<-x
    b<-c
  }
  if((h>=240)&(h<300)){
    r<-x
    b<-c
  }
  if((h>=300)&(h<360)){
    r<-c
    b<-x
  }

  m<-l-c/2
  return(clamp(c(r+m,g+m,b+m)))
}

3.5.4 Colorisation

Afin de coloriser une image, il suffit d’affecter à chaque pixel une teinte constante. Cette opération peut être réalisée en trois étapes :

  • Conversion RGB vers HSL
  • Modification de la composante H correspondant à la teinte désirée
  • Conversion HSL vers RGB

Exercice 8
Écrire une fonction R colorize(src,hue) qui colorise une image src passée en paramètre avec la teinte hue, correspondant à un angle exprimé en degrés.

display(colorize(img1, 350),colorize(img1, 200))

3.5.5 Gradient radial

Le principe du gradient radial consiste à diminuer la luminosité (lightness) d’un pixel en fonction de sa distance \(d\) au centre \((x_C,y_C)\) d’un cercle de rayon \(r\) : la nouvelle luminosité \(L'\) correspond au minimum entre l’ancienne lumonisité \(L\) du pixel et la valeur \(f(d)\).

\[ L' = min(L,f(d)) \] \[ \text{avec} \quad f(d) = max(0,1-d/r) \quad \text{et} \quad d = \sqrt{ (x - x_C)^2 + (y - y_C)^2 } \]

Comme pour la colorisation, cette opération peut être réalisée en 3 étapes pour chaque pixel :

  • Conversion RGB vers HSL
  • Calcul de la distance au centre du cercle
  • Calcul de la composante L selon la formule précédente
  • Conversion HSL vers RGB

Exercice 9
Écrire une fonction R gradient(src,centerx,centery,radius) qui applique le principe du gradient radial à une image src passée en paramètre pour un cercle de centre (centerx,centery) et de rayon radius.

display(img1,gradient(img1,90,60,52))

4 Transformations géométriques

Jusqu’à présent, nous avons utilisé les applications linéaires pour modifier les composantes d’un pixel. Avec les transformations géométriques, nous allons modifier les coordonnées \((x,y)\) d’un pixel afin de le déplacer au sein de l’image.

Toutes les transformations que nous allons implémenter se feront relativement à un centre passé en paramètre de chaque fonction. Par conséquent, avant d’appliquer la transformation, il faudra translater le pixel en lui retranchant les coordonnées du centre. Puis, après avoir appliqué la transformation, le pixel résultat sera à nouveau translaté dans le sens opposé en lui ajoutant les coordonnées du centre.

4.1 Changement d’échelle

Un changement d’échelle de centre \((x_C,y_C)\) et de rapport \(\lambda\) est défini comme suit :

\[ \left( \begin{array}{c} x' \\ y' \end{array} \right) = \begin{pmatrix} \lambda & 0 \\ 0 & \lambda \end{pmatrix} \left( \begin{array}{c} x - x_C \\ y - y_C \end{array} \right) + \left( \begin{array}{c} x_C \\ y_C \end{array} \right) \]

Exercice 10
Écrire une fonction R zoom(src,centerx,centery,factor) qui applique un changement d’échelle de rapport factor et de centre (centerx,centery) à une image src passée en paramètre.

display(img1,zoom(img1,width/2,height/2,0.5))

On constate qu’avec un facteur de changement d’échelle \(>1\), l’image résultat est creuse. En effet, en plus de déplacer le pixel, la fonction devrait affecter la même couleur aux pixels voisins dans un carré de coté factor.

display(img1,zoom(img1,width/4,height/4,4))

Exercice 10 bis
Écrire une fonction R zoomfill(src,centerx,centery,factor) qui corrige le défaut de la première fonction de changement d’échelle en procédant au remplissage des pixels voisins.

display(img1,zoomfill(img1,width/4,height/4,4))

Mais cette solution n’est pas entièrement satisfaisante. Pour éviter d’obtenir une image creuse, il suffit de parcourir les pixels de l’image destination et d’appliquer la transformation inverse pour retrouver le pixel antécédent. Cette dernière solution est la plus efficace.

Exercice 10 ter
Écrire une fonction R smartzoom(src,centerx,centery,factor) qui corrige le défaut de la deuxième fonction de changement d’échelle en appliquant la transformation inverse aux pixels de l’image destination.

display(img1,smartzoom(img1,width/4,height/4,4))

4.2 Rotation

Une rotation dans le plan de centre \((x_C,y_C)\) et d’angle \(\alpha\) est définie comme suit :

\[ \left( \begin{array}{c} x' \\ y' \end{array} \right) = \begin{pmatrix} \cos\alpha & -\sin\alpha \\ \sin\alpha & \cos\alpha \end{pmatrix} \left( \begin{array}{c} x - x_C \\ y - y_C \end{array} \right) + \left( \begin{array}{c} x_C \\ y_C \end{array} \right) \]

Exercice 11
Écrire une fonction R rotate(src,centerx,centery,angle) qui applique une rotation de centre (centerx,centery) et d’angle angle à une image src passée en paramètre.

display(img1,rotate(img1,width/2,height/2,pi/2))

Avec une autre valeur d’angle, on constate quelques jolis effets de Moiré. Ceci est dû à la nature discrète de la grille de pixels de l’image. Après transformation, les coordonnées du pixel sont modifiées pour n’en conserver que la partie entière, ce qui explique que certaines cases de l’image ne sont jamais remplies.

display(img1,rotate(img1,width/2,height/2,pi/3))

Exercice 11 bis
Écrire une fonction R smartrotate(src,centerx,centery,angle) qui corrige le défaut de la première fonction de rotation.

display(img1,smartrotate(img1,width/2,height/2,pi/3))

4.3 Twist

Le twist consiste en une rotation du plan de centre \((x_C,y_C)\), dont l’angle \(\alpha\) est fonction de la distance du pixel au centre.

\[ \left( \begin{array}{c} x' \\ y' \end{array} \right) = \begin{pmatrix} \cos\alpha & -\sin\alpha \\ \sin\alpha & \cos\alpha \end{pmatrix} \left( \begin{array}{c} x - x_C \\ y - y_C \end{array} \right) + \left( \begin{array}{c} x_C \\ y_C \end{array} \right) \]

\[\alpha = \rho \sqrt{ (x - x_C)^2 + (y - y_C)^2 }\]

Exercice 12
Écrire une fonction R twist(src,centerx,centery,rho) qui applique un twist de centre (centerx,centery) et de facteur rho à une image src passée en paramètre.

display(img1,twist(img1,width/2,height/2,0.01))

5 Récréation

5.1 Sauvegarde d’images

Le package PNG permet aussi d’enregistrer une image stockée en mémoire sous la forme d’un fichier au format PNG.

# Sauvegarde du tableau img dans le fichier image.png
writePNG(img,target="image.png");

5.2 Création d’un fichier vidéo à partir de plusieurs images

Pour finir, il est possible d’assembler plusieurs images PNG pour en faire un fichier vidéo. L’outil FFmpeg est particulièrement adapté pour ce genre de tâches : si vos images sont nommées de la façon suivante img001.png, img002.png, img003.png, … Vous pouvez taper la commande système ci-dessous dans un terminal.

Pour générer un fichier vidéo video.mp4 au format MP4 :

ffmpeg -framerate 25 -i img%03d.png -c:v libx264 -b:v 2M -vf format=yuv420p video.mp4

Pour générer un fichier vidéo video.webm au format WebM :

ffmpeg -framerate 25 -i img%03d.png -c:v libtheora -b:v 2M video.ogg

Exercice 13
Écrire une fonction R filename(count,path) qui retourne une chaine de caractères correspondant au nom de l’image numéro count dans le sous-répertoire path.

print(filename(42,"movie"))
## [1] "movie/img042.png"

Exercice 14
Écrire une fonction R rotozoom(src) qui génère une suite d’images PNG résultant de l’application d’une rotation suivie d’un zoom à l’image src passée en paramètre.

Exercice 15
Écrire une fonction R movie(src) qui génère une suite d’images PNG résultant de l’application de plusieurs filtres à l’image src passée en paramètre. Soyez créatifs !!!

6 Compression d’images

Le principe de la compression destructive consiste à réduire la taille d’un fichier image en acceptant de perdre de l’information et donc d’altérer l’image d’origine. La technique présentée par la suite et détaillée dans le document compression.pdf, permet de diviser la taille d’un fichier par 4.

La technique proposée sera uniquement appliquée à des images de luminance (i.e. en niveaux de gris). Pour des considérations de symétrie, la luminance d’un pixel sera représentée par un nombre \(l \in [-1,1]\) lors de l’étape de compression, \(-1\) correspondant au noir et \(1\) au blanc.

Une première version naïve de la compression consiste à diviser l’image en paquets carrés de 4 pixels, représentés par des vecteurs de \([-1,1]^4\). Au sein de chaque paquet, on identifie le pixel avec la plus grande luminance en valeur absolue. Seul ce pixel prépondérant est conservé lors de la compression, les luminances des autres pixels étant fixées à \(0\).

Pour supprimer l’effet dentelle engendré par \(3/4\) de pixels gris avec une luminance à \(0\), on peur dupliquer la luminance du pixel prépondérant dans les 3 autres.

Cette première version de la compression peut-être réalisée en 3 étapes :

Exercice 16
Écrire une fonction R compression(src) qui modifie une image de luminance src selon la méthode de compression naïve décrite précédemment.

imgg1<-greyscale(img1)
imggc1<-compression(imgg1)
display(imggc1, abs(imgg1-imggc1))

Afin d’atténuer l’effet de pixelisation, il paraît judicieux de privilégier une autre base de vecteurs dans \([-1,1]^4\) qui rend mieux compte des différences de luminance au sein d’un paquet de 4 pixels. Une version smart du même algorithme de compression consiste à :

Exercice 16 bis
Écrire une fonction R smartcompression(src) qui modifie une image de luminance src selon la méthode de compression smart décrite précédemment.

imgg1<-greyscale(img1)
imggsc1<-smartcompression(imgg1)
display(imggsc1, abs(imgg1-imggsc1))