Statistical analysis

The dimensions of the deformation are very high (1521 vectors in 3D, 4563 dimensions). In order to visualise the results in 2D or 3D, we can use a Principal Component Analysis (PCA). Here is an example based on the first three axes (interactive figure):

The code to plot the PCA:

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()