
Basic Procrustes Analysis


Igor
#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
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]

Forum

Support

Gallery
Igor Pro 9
Learn More
Igor XOP Toolkit
Learn More
Igor NIDAQ Tools MX
Learn More
Help would be appreciated.
My error. You are correct it compiles correctly
March 5, 2016 at 04:28 am - Permalink
A.G.
WaveMetrics, Inc.
March 4, 2016 at 10:28 am - Permalink