comparing shape of two 3D point clouds with Procrustes Analysis

geometry, r, shapes

Solution

Thanks to Julien Claude's book `Morphometrics with R` we have some handy code to do the same as that matlab function.

He provides some functions to calculate the full Procrustes distance which he defines as 'the square root of the sum of the squared distances between homologous coordinates of superimposed configurations (previously scaled to unit size)', just as the matlab function is defined.

# first, scale the coordinates to unit centroid size, and return both the scaled coords and the centroid size

centsiz<-function(M)
       {p<-dim(M)[1]
         size<-sqrt(sum(apply(M, 2,var))*(p-1))
         list("centroid_size" = size,"scaled" = M/size)}

# second, translate the coords so that its centroid is set at the origin

trans1<-function(M){scale(M,scale=F)}

# third, prepare the fPsup function to perform the full Procrustes superimposition of M1 onto M2. In the output, DF is the Full Procrustes distance between M1 and M2.

fPsup<-function(M1, M2) { 
       k<-ncol(M1)
          Z1<-trans1(centsiz(M1)[[2]])
          Z2<-trans1(centsiz(M2)[[2]])
          sv<-svd(t(Z2)%*%Z1)
          U<-sv$v; V<-sv$u; Delt<-sv$d
          sig<-sign(det(t(Z2)%*%Z1))
          Delt[k]<-sig*abs(Delt[k]) ; V[,k]<-sig * V[,k]
          Gam<-U%*%t(V)
          beta<-sum(Delt)
          list(Mp1=beta*Z1%*%Gam,Mp2=Z2,rotation=Gam,scale=beta,
                  DF=sqrt(1-beta^2))}

# test it out...
library(shapes) # so we can use the built-in data
data(gorf.dat) # Female gorilla skull data, 8 landmarks in 2 dimensions, 30 individuals

# calculate procrustes distance for individuals 1 and 2
fPsup(gorf.dat[,,1], gorf.dat[,,2])$DF
[1] 0.0643504

# Claude provides a check with a function that calculates the interlandmark distances between two configurations, which we can then sqrt the sum of to get the matlab-defined procrustes distance. 

ild2<-function(M1, M2){sqrt(apply((M1-M2)^2, 1, sum))}

# test it out...
test<-fPsup(gorf.dat[,,1], gorf.dat[,,2])
test$DF
[1] 0.0643504
sqrt(sum(ild2(test$Mp1, test$Mp2)^2))
[1] 0.0643504 # the same

If you just want to stick with the `shapes` package, the Riemannian shape distance function computes nearly identical results:

library(shapes)
riemdist(gorf.dat[,,1], gorf.dat[,,2])
[1] 0.0643949

UPDATE I've had some correspondence with Ian Dryden, the author of the `shapes` package. He wrote that to get the full Procrustes distance you just need to use `sin(riemdist)`. So the full Procrustes distance between the first two female gorillas is:

sin(riemdist(gorf.dat[,,1],gorf.dat[,,2])) 
[1] 0.0643504

And if we wanted to make our own function `fpdist` which does the same thing:

fpdist<-function(x, y, reflect = FALSE){
sin(riemdist(x,y,reflect=reflect))
}

fpdist(gorf.dat[,,1],gorf.dat[,,2]) 
[1] 0.0643504

Note that the gorilla data used above are 2D, but 3D data works fine also:

library(shapes) # so we can use the built-in data
data(macm.dat) # Male macaque skull data. 7 landmarks in 3 dimensions, 9 individuals

# calculate procrustes distance for macaque individuals 1 and 2
# Claude's method 1
fPsup(macm.dat[,,1], macm.dat[,,2])$DF
[1] 0.1215633

# Claude's method 2
test<-fPsup(macm.dat[,,1], macm.dat[,,2])
sqrt(sum(ild2(test$Mp1, test$Mp2)^2))
[1] 0.1215633

# using the shapes package
fpdist(macm.dat[,,1], macm.dat[,,2])
[1] 0.1215633

Is that what you were after?

Problem

I have two three-dimensional point clouds. I would like to compare their shapes and extent. I assume that Procrustes Analysis is the way to go. I have installed the package 'shapes' which provides several types of Procrustes Analysis, e.g. General Procrustes Analysis (GPA). I think, I am missing something here. What I was expecting a function to which I pass the two 3D matrices and it will return me a value on how good they match/correlate, e.g. a value between 0 - 1. Something like: ``` procrustes.distance(A,B) # A and B each being 3x100 ``` Basically something like procrustes in Matlab.

Original source