procOPA: Ordinary Procrustes analysis

Description

Ordinary Procustes analysis : the matching of one configuration to another using translation, rotation and (possibly) scale. Reflections can also be included if desired. The function matches configuration B onto A by least squares.

Usage

procOPA(A, B, scale = TRUE, reflect = FALSE)

Arguments

k x m matrix (or complex k-vector for 2D data), of k landmarks in m dimensions. This is the reference figure.

k x m matrix (or complex k-vector for 2D data). This is the figure which is to be transformed.

scale

logical indicating if scaling is required

reflect

logical indicating if reflection is allowed

Value

A list with components:
RThe estimated rotation matrix (may be an orthogonal matrix if reflection is allowed)
sThe estimated scale matrix
AhatThe centred configuration A
BhatThe Procrustes registered configuration B
OSSThe ordinary Procrustes sum of squares, which is $|Ahat-Bhat|^2$

References

Dryden, I.L. and Mardia, K.V. (1998). Statistical shape analysis. Wiley, Chichester.

Examples

Run this code

data(digit3.dat)

A<-digit3.dat[,,1]
B<-digit3.dat[,,2]
ans<-procOPA(A,B) 
plotshapes(A,B,joinline=1:13)
plotshapes(ans$Ahat,ans$Bhat,joinline=1:13)

#Sooty Mangabey data
data(sooty.dat)
A<-sooty.dat[,,1]   #juvenile
B<-sooty.dat[,,2]   #adult
par(pty="s")
plot(A,xlim=c(-2000,3000),ylim=c(-2000,3000))
lines(A[c(1:12,1),])
points(B)
lines(B[c(1:12,1),],lty=2)
#match B onto A
out<-procOPA(A,B)
#rotation angle
print(atan2(out$R[1,2],out$R[1,1])*180/pi)
#scale
print(out$s)
plot(A,xlim=c(-2000,3000),ylim=c(-2000,3000))
lines(A[c(1:12,1),])
points(out$Bhat)
lines(out$Bhat[c(1:12,1),],lty=2)
#match A onto B
out<-procOPA(B,A)
#rotation angle
print(atan2(out$R[1,2],out$R[1,1])*180/pi)
#scale
print(out$s)
plot(B,xlim=c(-2000,3000),ylim=c(-2000,3000))
lines(B[c(1:12,1),],lty=2)
points(out$Bhat)
lines(out$Bhat[c(1:12,1),])