1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
|
find_resolution<-function(canal_R,cote,method="mclust",affichage=TRUE,seuil=0.8)
{
##
# Input:
# ->cote: coté du carré en mm
# ->canal_R: Canal Rouge de l'image
# ->method: "kmeans" ou "hclust"
# ->seuil: pour seuiller l'image et prélever les 4 lasers
# Output:
# -> taille du pixel en mm
#points susceptibles d'appartenir au sommet du carré
pts=which(canal_R>max(canal_R)*seuil,arr.ind=TRUE)
ptm<-proc.time()
# Détermination des 4 groupes des nuages de points
N=4
if(method=="kmeans"){
cadre=matrix(c(1,1,dim(canal_R)[1],dim(canal_R)[1],1,dim(canal_R)[2],1,dim(canal_R)[2]),ncol=2)
cl<-tryCatch(kmeans(pts,cadre), error= function(e) {print(e);kmeans(pts,N)})
sommets=cl$centers
cl=cl$cluster
}else if(method=="hclust")
{
hca=hclust(dist(pts)^2)
cl=cutree(hca,k=N)
sommets=matrix(1:(N*2),ncol=2)*0
for(i in 1:N){
cluster=matrix(pts[cl==i],ncol=2)
sommets[i,]=colSums(cluster)/dim(cluster)[1]
}
}else if(method=="mclust"){
library(mclust)
cl=Mclust(pts)
nb_cl=cl$G
cl=cl$classification
nb_pts_par_cl=1:nb_cl*0 #allocation
for(i in 1:nb_cl){ nb_pts_par_cl[i]=sum(cl==i)}
ordre=sort(nb_pts_par_cl,decreasing=TRUE,index.return=TRUE) #classement par plus grand nombre de points
sommets=matrix(1:min(N*2,2*nb_cl),ncol=2)*0
for(i in 1:min(N,nb_cl)){
cluster=matrix(pts[cl==ordre$ix[i]],ncol=2)
sommets[i,]=colSums(cluster)/dim(cluster)[1] #prélèvement des 4 clusters
}
N=nb_cl
}else if(method=="minmax"){
N=0
cl1=c( min(pts[,1]),min(pts[,2]) )
cl2=c( max(pts[,1]),max(pts[,2]) )
sommets=rbind(cl1,cl2)
}else{
print("Mauvaise entrée pour l'argument: method (hclust,kmeans,mclust,minmax) ")
return (0)
}
#Détermination d'une diagonale <=> distance la plus grande des clusters
ind_clust_diag=as.matrix(dist(sommets))
clust_diag= arrayInd(which.max(ind_clust_diag),dim(ind_clust_diag))
centre_carre= colSums(sommets[clust_diag,]) /2 #milieux de la diagonale
distance_moyenne=mean( sqrt(rowSums(( cbind(sommets[,1]-centre_carre[1],sommets[,2]-centre_carre[2] ) )^2)) )
cote_pixel=distance_moyenne*2/sqrt(2)
taille_pixel=cote/cote_pixel
print(proc.time()-ptm )
# Affichage
if(affichage==TRUE)
{
windows()
layout(matrix(1:4,2,2))
image(canal_R,col=grey(0:255/255))
title(c("Image faveur Rouge:"," 2/3.Red-1/3(Green+Blue)"))
image(canal_R>max(canal_R)*seuil,col=grey(0:255/255))
title("image seuillé")
plot(pts[,1],pts[,2])
if(N>0){for(i in 1:N){points(matrix(pts[cl==i],ncol=2),col=i) } }
points (sommets[,1],sommets[,2],pch=3,col=N+1)
title(paste("Clustering par la méthode:",method))
mtext(paste("taille pixel=",taille_pixel," mm"))
}
return(taille_pixel)
} |
Partager