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]