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
		print "theta=",theta, "rad"
	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

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]

procrustes(wave1,wave2)

Display w1[][1] vs w1[][0]
AppendToGraph w2[][1] vs w2[][0]