# Basic Procrustes Analysis

#pragma rtGlobals=3     // Use modern global access method and strict wave access.

// Very basic example of Procrustes Analysis.
// wave1 and wave2 are expected to have two columns representing
// X and Y components.
Function procrustes(wave1,wave2)
Wave wave1,wave2

Variable rows1=dimsize(wave1,0)
Variable cols1=dimsize(wave1,1)
Variable rows2=dimsize(wave2,0)
Variable cols2=dimsize(wave2,1)

// Step 1:
// do not modify the original waves.  Create local waves w1
// and w2.  Subtract the mean from each column.  That shifts
// each shape to its origin.
MatrixOP/O w1=subtractMean(wave1,1)
MatrixOP/O w2=subtractMean(wave2,1)

// Step 2: perform scaling by dividing by stdv.
//  Get the scaling factors and arrange as col wave.
MatrixOP/O/FREE stdv1=sqrt(varcols(w1))^t
MatrixOP/O/FREE stdv2=sqrt(varcols(w2))^t
//  Scale
MatrixOP/O w1=scaleCols(w1,rec(stdv1))
MatrixOP/O w2=scaleCols(w2,rec(stdv2))

// test here:
MatrixOP/O/FREE v1=sum(varCols(w1)+varCols(w2))
if(v1[0]!=4)
Print "Bad scaling factors."
endif

// At this point we removed position and scale.  Next comes rotation.
if(rows1==rows2)
MatrixOP/O/FREE x1=col(w1,0)
MatrixOP/O/FREE x2=col(w2,0)
MatrixOP/O/FREE y1=col(w1,1)
MatrixOP/O/FREE y2=col(w2,1)
MatrixOP/O yc=sum(x2*y1-y2*x1)
MatrixOP/O xc=sum(x2*x1+y2*y1)
Variable theta=atan2(yc[0],xc[0])
// create the rotation matrix:
Make/O/D/N=(2,2) rotMat
set2DRotMatrix(rotMat,theta)

// find the initial procrustes distance
variable dist=getProcrustesDist(w1,w2)
print " before dist=",dist

// Step 3. rotate w2 to match w1:
MatrixOP/O w2=(rotMat x w2^t)^t
dist=getProcrustesDist(w1,w2)
print " after dist=",dist
endif
End

Function set2DRotMatrix(rotMat,theta)
Wave rotMat
Variable theta

rotMat[0][0]=cos(theta)
rotMat[1][1]=cos(theta)
rotMat[0][1]=-sin(theta)
rotMat[1][0]=sin(theta)
End

Function getProcrustesDist(w1,w2)
Wave w1,w2

Matrixop/o/free aa=sum(powr(w1-w2,2))
return sqrt(aa[0])
End

Here is a simple example inspired by the "big W":
Make/N=(5,2) wave1,wave2
wave1[0][0]= {0,0.326797,0.0588235,0.326797,0}
wave1[0][1]= {1,0.703947,0.486842,0.217105,0}
wave2[0][0]= {1.9281,1.59477,1.86275,1.62092,1.87582}
wave2[0][1]= {0.263158,0.375,0.559211,0.756579,0.901316}
Display wave1[][1] vs wave1[][0]
AppendToGraph wave2[][1] vs wave2[][0]

Now execute the function above:
procrustes(wave1,wave2)

Display the results:
Display w1[][1] vs w1[][0]
AppendToGraph w2[][1] vs w2[][0]

I couldn't get this to compile. The problem line seems to be: set2DRotMatrix(rotMat,theta)
Help would be appreciated.

My error. You are correct it compiles correctly

Forum

Support

Gallery