r - 使用 Procrustes Analysis 比较两个 3D 点云的形状

Question

我有两个三维点云。我想比较一下它们的形状和范围。我假设 Procrustes Analysis是要走的路。我已经安装了包'shapes ' 它提供了几种类型的 Procrustes 分析，例如一般 Procrustes 分析 (GPA)。我想，我在这里遗漏了一些东西。我期待一个函数，我将两个 3D 矩阵传递给它，它会返回一个关于它们匹配/相关程度的值，例如一个介于 0 - 1 之间的值。类似于:

procrustes.distance(A,B) # A and B each being 3x100 

基本上类似于 procrustes在Matlab中。

Answer 1

多亏了 Julien Claude 的书Morphometrics with R，我们有了一些方便的代码来执行与该 matlab 函数相同的操作。

他提供了一些函数来计算完整的 Procrustes 距离，他将其定义为“叠加配置的同源坐标之间的平方距离之和的平方根(之前缩放为单位大小)”，就像 matlab 函数的定义一样.

# 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

如果您只想坚持使用 shapes 包，黎曼形状距离函数计算出几乎相同的结果:

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

更新 我与 shapes 包的作者 Ian Dryden 有一些通信往来。他写道，要获得完整的 Procrustes 距离，您只需要使用 sin(riemdist)。所以前两只雌性 gorilla 之间的完整 Procrustes 距离是:

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

如果我们想创建自己的函数 fpdist 来做同样的事情:

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

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

请注意，上面使用的 gorilla 数据是 2D，但 3D 数据也可以正常工作:

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

这就是你想要的吗？