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