Traitement et analyse d'images
en langage R
L'essentiel de cette page
Voici une nouvelle page qui se fabrique. R n'est absolument pas fait pour le traitement et l'analyse d'images (TAI) mais ça peut être un gain de temps certains si on combine analyses statistiques et TAI dans le même environnement (R) sans passer par plusieurs logiciels.
1) Comment visualiser une image ?
Il existe de nombreuses fonctions pour voir les images (TIF, JPG, PNG...).
Prenons l'exemple d'une image png.
library("png")
img1<-readPNG("png.png")
Pour la visualiser, on pourra faire appelle à la fonction (très pratique) display() développée par Frédéric Blanchard.
display(img1) # cf. ci-après pour disposer de la fonction display
Code à copier-coller pour bénéficier de la fonction display().
display<-function(...){
imageList<-list(...) ; totalWidth<-0 ; maxHeight<-0 ;
for (img in imageList){ if(is.character(img)) img<-readPNG(img) ; dimg<-dim(img) ; totalWidth<-totalWidth+dimg[2] ; maxHeight<-max(maxHeight,dimg[1]) }
par(mar=c(0,0,0,0)) ; plot(c(0,totalWidth),c(0,maxHeight),type="n",asp=1,xaxt="n",yaxt="n",xlab="x",ylab="y") ; offset<-0 ;
for (img in imageList){ dimg<-dim(img) ; rasterImage(img,offset,0,offset+dimg[2],dimg[1])
offset<-offset+dimg[2] }}
2) Visualiser les caractéristiques de l'image (répartition des couleurs par intensité)
Un histogramme permet de se donner une idée de la répartition de l'intensité des valeurs sauf qu'il faut pouvoir séparer les composantes rouge, verte et bleue (rgp).
Voici un code qui permet d'avoir accès à une fonction d'histogramme.
hist_img(img1)
Filtrer une image pour n'afficher en couleur ou en niveau de gris qu'une ou plusieurs composantes parmi les composantes RVB (Rouge, Vert, Bleu)
La fonction filter_img() permet d'extraire une composante (rouge, vert ou bleu) d'une image ou plusieurs (rouge et vert à la fois par ex). Cette composante peut être convertie automatiquement en niveau de gris.
Rouge
img2 <- filter_img(img1,type="r") ; display(img2)
Vert
img2 <- filter_img(img1,type="v") ; display(img2)
Bleu
img2 <- filter_img(img1,type="b") ; display(img2)
filter_img <- function(img,type="r",grey=F,compare=T,return=T){
if (grey==F) {
temp <- img
if (type == "r"){temp [,,2] <- 0 ; temp[,,3] <- 0
} else if ((type == "v")|(type == "g")){temp [,,1] <- 0 ; temp[,,3] <- 0
} else if (type == "b"){temp [,,1] <- 0 ; temp[,,2] <- 0
} else if ((type == "rv")|(type == "vr")|(type == "gr")|(type == "rg")){temp[,,3] <- 0
} else if ((type == "rb")|(type == "br")){temp[,,2] <- 0
} else if ((type == "vb")|(type == "bv")|(type == "gb")|(type == "bg")){temp[,,1] <- 0 }
} else if (grey==T) {
if (type == "r"){temp <- img[,,1]
} else if ((type == "v")|(type == "g")){temp <- img[,,2]
} else if (type == "b"){temp <- img[,,3]
} else if ((type == "rv")|(type == "vr")|(type == "rg")|(type == "gr")){temp<- (img[,,1]+img[,,2])/2
} else if ((type == "rb")|(type == "br")){temp<- (img[,,1]+img[,,3])/2
} else if ((type == "vb")|(type == "bv")|(type == "gb")|(type == "bg")){temp<- (img[,,2]+img[,,3])/2}}
if (return == T) {return(temp)}}
2.2) Extraire le niveau de gris
###############################
# NVEAU DE GRIS
###############################
gray_scale_img <- function(img){
temp_img <- img[,,1]*0.2126+img[,,2]*0.7152+img[,,3]*0.0722
return(temp_img)
}
img4 <- gray_scale_img(img1)
display(img1,img4)
2.3) Afficher une composante sans filtrer (rapide)
background <- function(img, type="r") {
temp <- img
if (type=="r") {temp[,,2] <- 0 ; temp[,,3] <- 0
}else if(type=="g") {temp[,,1] <- 0 ; temp[,,3] <- 0
}else if(type=="b") {temp[,,1] <- 0 ; temp[,,2] <- 0}
display(temp)
return(temp)
}
background(img1,type="r")
2.4) Normaliser les composantes de l'image (expansion dynamique)
Objectif : forcer la valeur min à être à 0 et la valeur max à 255.
normalize_img <- function(img) {
temp_img <- img
temp_img[,,1] <- temp_img[,,1] - min(temp_img[,,1])
temp_img[,,1] <- temp_img[,,1] / max(temp_img[,,1])
temp_img[,,2] <- temp_img[,,2] - min(temp_img[,,2])
temp_img[,,2] <- temp_img[,,2] / max(temp_img[,,2])
temp_img[,,3] <- temp_img[,,3] - min(temp_img[,,3])
temp_img[,,3] <- temp_img[,,3] / max(temp_img[,,3])
return(temp_img)
}
img4 <- normalize_img(img1)
Avant
Après (contraste augmenté sans saturation du signal)
3) Réduire le nombre de pixels d'une image afin d'accélérer les calculs et de définir le bon procédé de traitement avant de revenir à une image pleine résolution.
En projet : extraire une matrice moyenne réduite. Bientôt disponible sur ce site avec une fonction de binarisation. Un filtre moyen. Un filtre médian. Un filtre de convolution. Une fonction d'érosion. A bientôt.
3.1) Réduire le nombre de pixels d'une image pour appliquer des traitements rapides
La fonction redim_img() a été programmée pour réduire la taille d'une image. Il suffit d'appliquer un coefficient (coef).
On peut aussi forcer la taille de l'image avec height ou width.
Enfin, un lissage est prévu, mais non encore ajouté avec les paramètre smooth et s.type qui devront être programmés dans une version ultérieure.
redim_img <- function(img,coef=c(),width=c(),height=c(),smooth=F,s.type = mean) {
if(length(coef)>0) {
if (coef[1] == 1){cat("Le coefficient de redimension doit être inférieur ou supérieur à 1.\n")
} else if (coef[1] < 1) {
# Calcul des dimensions de la nouvelle image avec conservation maximale des proportions
height_img <- nrow(img)
width_img <- ncol(img)
color <- dim(img)[3]
prop <- width_img/height_img
height_new <- coef[1]*height_img
width_new <- coef[1]*width_img
diff1 <- ceiling(width_new)/ceiling(height_new)-prop
diff2 <- floor(width_new)/floor(height_new)-prop
diff3 <- ceiling(width_new)/floor(height_new)-prop
diff4 <- floor(width_new)/ceiling(height_new)-prop
diff <- abs(c(diff1,diff2,diff3,diff4))
ind <- which(diff==min(diff))
if (ind[1]==1) {height_new <- ceiling(height_new); width_new <- ceiling(width_new)
} else if (ind[1]==2) {height_new <- floor(height_new); width_new <- floor(width_new)
} else if (ind[1]==3) {height_new <- floor(height_new); width_new <- ceiling(width_new)
} else if (ind[1]==4) {height_new <- ceiling(height_new); width_new <- floor(width_new)
}
if (smooth == F) {
wd <- round(seq(1,width_img,length.out=width_new),0)
ht <- round(seq(1,height_img,length.out=height_new),0)
ind <- expand.grid(ht,wd)
temp <- array(0,dim=c(height_new,width_new,color)) ; k= 0
for (i in 1:width_new) {
for (j in 1:height_new) {
k <- k+1
temp[j,i,]<- img[ind[[1]][k],ind[[2]][k],1:color]
}
}
return(temp)
}
}else if (coef[1] > 1) {
}
} else if (length(width)>0){
}
}
img2 <- redim_img(img1,coef=0.1)
display(img2)
3.2) Filtre médian et filtre moyen
#############################
# GENERER UN FOND MEDIAN à 5 points
#############################
median_img <- function(img, taille=40) {
temp_img <- img
largeur <- ncol(img) ; hauteur <- nrow(img)
pas <- round(largeur/100)+1
compteur =c(1:100)*pas
i = 1
compteur =c(1:80)*25
i = 1
for (nc in c(1:largeur)) {
if (nc == compteur[i]) {
percent = round(compteur[i]/pas)
cat(percent,"%\n")
i = i+1
}
for (nl in c(1:hauteur)) {
#medium_img[nl,nc,1] <- median(img1[nl-pas:nl+pas,nc-pas:nc+pas,1] )
#medium_img[nl,nc,2] <- median(img1[nl-pas:nl+pas,nc-pas:nc+pas,2] )
#medium_img[nl,nc,3] <- median(img1[nl-pas:nl+pas,nc-pas:nc+pas,3]
if (((nl-taille)<0)|((nl+taille)>hauteur)|((nc-taille)<0)|((nc+taille)>largeur)){
itt <- 0
} else {itt <- taille}
temp_img[nl,nc,1] <- median(c(img[c(nl-itt,nl+itt),c(nc-itt,nc+itt),1],img[nl,nc,1] ))
temp_img[nl,nc,2] <- median(c(img[c(nl-itt,nl+itt),c(nc-itt,nc+itt),2],img[nl,nc,2] ))
temp_img[nl,nc,3] <- median(c(img[c(nl-itt,nl+itt),c(nc-itt,nc+itt),3],img[nl,nc,3]))
}}
return(temp_img)
}
medium_img <- median_img(img1)
x11()
display(img1,medium_img)
###############################
# AVERAGE
###############################
filtre_average <- function(img,taille=3) {
temp_img <- img
largeur <- ncol(img)
hauteur <- nrow(img)
pas <- round(largeur/100)+1
compteur =c(1:100)*pas
i = 1
for (nc in c(1:largeur)) {
if (nc == compteur[i]) {
percent = round(compteur[i]/pas)
cat(percent,"%\n")
i = i+1
}
for (nl in c(1:hauteur)) {
if (((nl-taille)<0)|((nl+taille)>hauteur)|((nc-taille)<0)|((nc+taille)>largeur)){
itt <- 0
} else {itt <- taille}
temp_img[nl,nc,1] <- mean(img[c(nl-itt,nl+itt),c(nc-itt,nc+itt),1] )
temp_img[nl,nc,2] <- mean(img[c(nl-itt,nl+itt),c(nc-itt,nc+itt),2] )
temp_img[nl,nc,3] <- mean(img[c(nl-itt,nl+itt),c(nc-itt,nc+itt),3] )
}}
return(temp_img)}
red <- filtre_average(background(img2,type="r"))
display(red)
3.3) Binarisation d'image
binary <- function(img,seuil_rvb=c(0,0,0)) {
seuil = 0.4
temp <- img
seuil=0.4
temp[,,1][temp[,,1]>seuil] <- 1
temp[,,1][temp[,,1]<seuil] <- 0
temp[,,2][temp[,,2]>seuil] <- 1
temp[,,2][temp[,,2]<seuil] <- 0
temp[,,3][temp[,,3]>seuil] <- 1
temp[,,3][temp[,,3]<seuil] <- 0
temp2 <- temp[,,1]*1/3+temp[,,2]*1/3+temp[,,3]*1/3
temp2[temp2>0.3] <- 1
temp2[temp2<0.3] <- 0
display(temp2)
return(temp2)
}
binary(img2) -> imgb
display(imgb)
On remarquera que cette fonciton n'est pas encore complètement aboutie (manque de temps) : seuil (valeur fixée actuelement à 0,4) devrait être remplacé par seuil_rvb afin de pouvoir paramétrer manuellement le seuil de binarisation (idéalement avec un outil interactif).
Mots clefs : binarisation, threshold
3.4) Filtre d'érosion et de convolution
Comme ces fonctions ne sont pas toutes rapides et aboutis, voici 4 versions à tester (manque de temps pour les parachever).
#######################
# VERSION 01
#######################
erode <- function(img,win=5,neig_min=0.75,forme=0){
debut <- ceiling(win/2)
pas <- (win-debut)
largeur <- ncol(img)
longueur <- nrow(img)
img_temp <- img
for (col in c(3:(largeur-pas))){
for (row in c(3:(longueur-pas))){
points_periph <- img[(row-pas):(row+pas),(col-pas):(col+pas)]
points_periph <- points_periph[-debut,-debut]
points_periph <- as.vector(points_periph)
points_periph <- sort(points_periph)
effectif <- table(points_periph)/length(points_periph)
indices <- as.numeric(names(effectif))
freq <- effectif[indices==forme]
freq <- ifelse(length(freq)<1,0,freq)
img_temp[row,col] <- ifelse(freq < neig_min,indices[indices!=forme],forme)
}
}
return(img_temp)
}
#######################
# VERSION 02
#######################
erode <- function(img,neig_min=0.75,win=5,forme=0){
forme_oppose <- ifelse(forme==0,1,0)
debut <- ceiling(win/2)
pas <- win-debut
largeur <- ncol(img)
longueur <- nrow(img)
img_temp <- img
pixels <- win^2-1
indices <- which(imgb==forme,arr.ind=T)
ind <- which(indices[,1]>2&indices[,1]<(longueur-pas)&indices[,2]>2&indices[,2]<(largeur-pas))
indices <- data.frame(row=indices[ind,1],col=indices[ind,2])
nindices <- nrow(indices)
for (i in (1:nindices)) {
lesrow <- seq(unlist(indices[i,][1]-pas),unlist(indices[i,][1]+pas),1)
lescol<-seq(unlist(indices[i,][2]-pas),unlist(indices[i,][2]+pas),1)
points_periph <- img[lesrow,lescol]
points_periph <- as.vector(points_periph[-debut,-debut])
if(length(points_periph[points_periph==forme])>0){
effectif <- length(points_periph[points_periph==forme])/pixels
if (effectif < neig_min) {
img_temp[unlist(indices[i,][1]),unlist(indices[i,][2])]<-forme_oppose
}
} else {
img_temp[unlist(indices[i,][1]),unlist(indices[i,][2])]<-forme_oppose
}
temp <- c(temp,effectif)
}
return(img_temp)
}
imgc <- erode(imgb,0.3)
#######################
# VERSION 03 - A FINIR - FENETRE REGLABLE...
#######################
erode <- function(img,neig_min=0.75,win=5,forme=0){
# Base commune
largeur <- ncol(img)
longueur <- nrow(img)
pixels <- win^2
centre <- ceiling(win/2)
pas <- win-centre
img_temp <- img
if (forme == 0) {
forme_oppose <- 1
indice <- round(pixels*(1-neig_min),0)
print(indice)
indices <- which(img==forme,arr.ind=T)
ind <- which(indices[,1]>2&indices[,1]<(longueur-pas)&indices[,2]>2&indices[,2]<(largeur-pas))
ind_x <- unlist(indices[ind,1])
ind_y <- unlist(indices[ind,2])
n_ind <- length(ind)
for (i in (1:n_ind)) {
lesrow <- seq((ind_x[i]-pas),(ind_x[i]+pas),1)
lescol<-seq((ind_y[i]-pas),(ind_y[i]+pas),1)
somme <- sum(img[lesrow,lescol])
if (somme > indice) {
img_temp[ind_x[i],ind_y[i]] <- 1}
}
} else if (forme == 1) {
}
return(img_temp)
}
#######################
# VERSION 04 - A FINIR - fenetre 3*3 obligatoire - NE fonctionne pas.
#######################
erode <- function(img,neig_min=0.75,win=3,forme=0){
# Base commune
largeur <- ncol(img)
longueur <- nrow(img)
pixels <- win^2
centre <- ceiling(win/2)
pas <- win-centre
img_temp <- img
if (forme == 0) {
forme_oppose <- 1
indice <- round(pixels*(1-neig_min),0)
print(indice)
indices <- which(img==forme,arr.ind=T)
ind <- which(indices[,1]>2&indices[,1]<(longueur-pas)&indices[,2]>2&indices[,2]<(largeur-pas))
ind_x <- unlist(indices[ind,1])
print(ind_x)
ind_y <- unlist(indices[ind,2])
n_ind <- length(ind)
print("error1")
l_1 <- img[(ind_x-1),(ind_y-1)]
print(l_1)
l_2 <- img[ind_x,ind_y-1]
l_3 <- img[ind_x+1,ind_y-1]
l_4 <- img[ind_x-1,ind_y]
l_5 <- img[ind_x,ind_y]
l_6 <- img[ind_x+1,ind_y]
l_7 <- img[ind_x-1,ind_y+1]
l_8 <- img[ind_x,ind_y+1]
l_9 <- img[ind_x+1,ind_y+1]
print("error2")
temp_tab <- data.frame(l_1,l_2,l_3,l_4,l_5,l_6,l_7,l_8,l_9)
temp_result <- apply(temp_tab,1,sum)
temp_result <- ifelse(temp_result>indice,1,0)
img_temp[ind_x,ind_y]<- temp_result
return(img_temp)
} else if (forme == 1) {
}
}
imgc <- erode(imgb,0.9)
filtre_convolution <- function(img,taille=3) {
filtre <- c(-2,-1,0,-1,1,1,0,1,2)
filtre <- c(0,-1,0,-1,5,-1,0,-1,0)
# matrice de 2 lignes sur 3 colonnes
filtre <- matrix(filtre,nrow=3,ncol=3)
#decallage = (taille-1)/2+1
temp_img <- img
largeur <- ncol(img)
hauteur <- nrow(img)
pas <- round(largeur/100)+1
compteur =c(1:100)*pas
i = 1
for (nc in c(2:(largeur-1))) {
if (nc == compteur[i]) {
percent = round(compteur[i]/pas)
cat(percent,"%\n")
i = i+1
}
for (nl in c(2:(hauteur-1))) {
temp_img[nl,nc,1] <- sum(img[(nl-1):(nl+1),(nc-1):(nc+1),1]*filtre)
temp_img[nl,nc,2] <- sum(img[(nl-1):(nl+1),(nc-1):(nc+1),2]*filtre)
temp_img[nl,nc,3] <- sum(img[(nl-1):(nl+1),(nc-1):(nc+1),3]*filtre)
}
temp_img[temp_img<0]<-0
temp_img[temp_img>1]<-1
}
return(temp_img)
}
4) Identifier des catégories de pixels par ACP et k-means
df <- data.frame(R = as.vector(img1[,,1]),
V = as.vector(img1[,,2]),
B = as.vector(img1[,,3]))
pca <- prcomp(df)
pca$x[,1]<-pca$x[,1]-min(pca$x[,1])
pca$x[,1]<-pca$x[,1]*1/max(pca$x[,1])
hist(pca$x[,1])
img2 <- img1
img2[,,1] <- array(as.vector(t(pca$x[,1])), dim = c(nrow(img1), ncol(img1), 1))
img2[,,2] <- array(as.vector(t(pca$x[,1])), dim = c(nrow(img1), ncol(img1), 1))
img2[,,3] <- array(as.vector(t(pca$x[,1])), dim = c(nrow(img1), ncol(img1), 1))
display(img2)
Toutes les composantes (R, G et B ont été réduit à une seule par ACP (en récupérant ainsi un maximum de signal propreà l'image) : on peut ainsi ajouter d'autres signaux tel que du IR.
# Soit img1 une image png au format array à 3 dimensions
# NOMBRE DE CATEGORIES
X <- 5
df <- data.frame(R = as.vector(img1[,,1]),
V = as.vector(img1[,,2]),
B = as.vector(img1[,,3]))
centres_x = sample(df$R,X); centres_y = sample(df$V,X);centres_z = sample(df$B,X)
centers = data.frame(centres_x, centres_y, centres_z)
# Regrouper les points autour de ces centres, recalculer de nouveaux centres : itérations
cluster_1 <- kmeans(df , centers, 10) # 5 représente le nombre maxima d'iterations
summary(cluster_1) # pour obtenir une description de l'objet ainsi créé
col <- cluster_1$cluster-min(cluster_1$cluster)
col <- col*1/max(col)
img2 <- img1
#img2 <- apply(img1, 3, function(x) as.numeric(factor(x, levels = cluster_1$centers[,x]))).
type <- factor(cluster_1$cluster)
levels(type) <- cluster_1$centers[,1] ; type <- as.numeric(as.character(type))
img2[,,1] <- array(as.vector(t(type)), dim = c(nrow(img2), ncol(img2), 1))
type <- factor(cluster_1$cluster)
levels(type) <- cluster_1$centers[,2]; type <-as.numeric(as.character(type))
img2[,,2] <- array(as.vector(t(type)), dim = c(nrow(img1), ncol(img1), 1))
type <- factor(cluster_1$cluster)
levels(type) <- cluster_1$centers[,3]; type <- as.numeric(as.character(type))
img2[,,3] <- array(as.vector(t(type)), dim = c(nrow(img1), ncol(img1), 1))
img2 <- array(as.vector(t(as.numeric(img2))), dim = c(nrow(img2), ncol(img2), 3))
display(img2)
Ainsi une catégorisation automatique a été faite (ici 5 couleurs de pixels, on peut aussi réduire ou augmenter le nombre de catégories et intégrer des données tel que du IR).
Quand X = 2