/* This program computes true Procrustes distances from an n by p*k matrix of
specimens, where n is the number of specimens and p*k is the number of coordinates 
(i.e., p=landmarks * k=dimensions). Coordinates of the first specimen are in ROW 1, 
coordinates of the second are in ROW 2, etc. The output matrix reports Procrustes 
distance as the angle (rho) in radians between each pair of specimens. */

proc iml;
use two; read all into X; read all var _CHAR_ into name;
k=3; n=nrow(X); p=(ncol(X)/k); 
dmat=j(n,n);
y=j(1,ncol(X));

do f=1 to n; do j=1 to n;
  X1=X[f,]; X1=shape(X1,p,k);
  X2=X[j,]; X2=shape(X2,p,k);

  I=I(p); CaP=j(p,p,1/p);
  s1=sqrt(trace((I-CaP)*X1*X1`*(I-CaP)));
  s2=sqrt(trace((I-CaP)*X2*X2`*(I-CaP)));
  X1prime=(I-CaP)*X1/s1; X2prime=(I-CaP)*X2/s2;
  cov2=X1prime`*X2prime;
  CALL SVD (U,sig,V,cov2);
  sig2=abs(sig); S=sig/sig2; S=S`; S=diag(S);
  H=V*S*U`;
  X2star=X2prime*H;

  y=X1prime-X2star;
  d=sqrt(trace(y`*y));
  rho=2*arsin(d/2);
  dmat[j,f]=rho;
  if j=f then do; dmat[j,f]=0; end;
end; end;

print dmat[rowname=name colname=name];
create procd from dmat; append from dmat;
abort;