Analyse des déformations

Les déformations obtenues sont de dimensions importantes (dans cet exemple, 1521 vecteurs 3D, soit de dimensions 4563). Afin de les visualiser dans une figure en 2D ou 3D, nous pouvons utiliser une analyse en composantes principales (ACP). Voici un graphique de la répartition des points en traçant les 3 premières composantes (figure interactive) :

Voici le code permettant de calculer l’ACP :

library(ade4)
library(rgl)
library(geometry)

Nsu = 12 # number of subjects (pooled)
sigmaV = 0.1 # standard deviation of the deformation kernel

#open moments (MOM_final.txt)
aux=as.matrix(read.table('../MOM_final.txt',header=F,sep=" "))
NBM = aux[1,]
aux = aux[-1,]
MOM = matrix(0,Nsu,3*NBM[2])
for (k in 1:Nsu){
  a = (k-1)*NBM[2]+1
  b = (k-1)*NBM[2]+NBM[2]
  MOM[k,] = c(t(aux[a:b,]))}

#open control points (CP_final.txt)
CP=t(as.matrix(read.table('../CP_final.txt',header=F,sep=" ")))
CPx = CP[rep(1,NBM[2]),]
CPy = CP[rep(2,NBM[2]),]
CPz = CP[rep(3,NBM[2]),]
K = exp(-( (CPx - t(CPx))^2 + (CPy - t(CPy))^2 + (CPz - t(CPz))^2 )/sigmaV^2);
V1 = K %*% t(MOM[,seq(1, dim(MOM)[2], by=3)])
V2 = K %*% t(MOM[,seq(2, dim(MOM)[2], by=3)])
V3 = K %*% t(MOM[,seq(3, dim(MOM)[2], by=3)])
V = cbind(t(V1),t(V2),t(V3))

label = c(rep(2,times=4),rep(3,times=4),rep(4,times=4))

library(ade4)

gp=dudi.pca(V,scann=F,nf=10000,center=T,scale=F)
score = gp$li

plot3d(score[,1],score[,2],score[,3],size=10,aspect=1,box=0,axes=FALSE,col=label,xlim=range(min(score[,1]),max(score[,1])),ylim=range(min(score[,2]),max(score[,2])),zlim=range(min(score[,3]),max(score[,3])),xlab="PC 1",ylab="PC 2",zlab="PC 3")

bg3d(color="white")
par3d(FOV=0)
aspect3d("iso")

rgl.lines(c(min(score[,1]), max(score[,1])), c(min(score[,2]), min(score[,2])), c(min(score[,3]), min(score[,3])), color = "black")
rgl.lines(c(min(score[,1]), min(score[,1])), c(min(score[,2]), max(score[,2])), c(min(score[,3]), min(score[,3])), color = "black")
rgl.lines(c(min(score[,1]), min(score[,1])), c(min(score[,2]), min(score[,2])), c(min(score[,3]), max(score[,3])), color = "black")

couleur = palette()[(2)]
couleuralpha=rgb(col2rgb(couleur)[1,]/255,col2rgb(couleur)[2,]/255,col2rgb(couleur)[3,]/255,0.4)
splitModerns = split(score[,c(1,2,3)], as.factor(label))
appliedModerns = t(convhulln(splitModerns[[1]]))
pts = as.matrix(splitModerns[[1]])
rgl.triangles(pts[appliedModerns,1],pts[appliedModerns,2],pts[appliedModerns,3],col=couleuralpha,alpha=.2)

couleur = palette()[(3)]
couleuralpha=rgb(col2rgb(couleur)[1,]/255,col2rgb(couleur)[2,]/255,col2rgb(couleur)[3,]/255,0.4)
splitModerns = split(score[,c(1,2,3)], as.factor(label))
appliedModerns = t(convhulln(splitModerns[[2]]))
pts = as.matrix(splitModerns[[2]])
rgl.triangles(pts[appliedModerns,1],pts[appliedModerns,2],pts[appliedModerns,3],col=couleuralpha,alpha=.2)

couleur = palette()[(4)]
couleuralpha=rgb(col2rgb(couleur)[1,]/255,col2rgb(couleur)[2,]/255,col2rgb(couleur)[3,]/255,0.4)
splitModerns = split(score[,c(1,2,3)], as.factor(label))
appliedModerns = t(convhulln(splitModerns[[3]]))
pts = as.matrix(splitModerns[[3]])
rgl.triangles(pts[appliedModerns,1],pts[appliedModerns,2],pts[appliedModerns,3],col=couleuralpha,alpha=.2)
rglwidget()