/* PROGRAM: eff.sas * Kathy Gerber, University of Virginia * * Compute and display outline approximation based on EFF. * Original Fortran source code developed by F. James Rohlf * and Scott Fersonm Department of Ecology and Evolution, * State University of New York, Stony Brook, NY, 1992. * * These methods are based on: * 1. Ferson, S. F., F. J. Rohlf, and R. K. Koehn. 1985. * Measuring shape variation of two-dimensional outlines. * Syst. Zoology, 34:59-68. * * 2. Kuhl, F P and C R Giardina 1982. Elliptic fourier * features of a closed contour. Computer graphics and image * processing, 18:236-258. * * ASSUMES: Data has the form id x1 y1 x2 y2... xn yn * Pairs is the number of xi yi pairs, * Harmonics is the number of harmonics < pairs/2 * */ %LET pairs = 256 ; %LET harmonics = 64 ; %LET datasrc = "A Two" ; FILENAME myfile 'g:\shapes\efa\ktest.dta' lrecl=2103; %LET dsn = mytest ; %MACRO readdata ; %LET ptnames = ; %DO i = 1 %TO &pairs ; %LET ptnames= &ptnames x&i y&i ; %END ; ; RUN ; DATA &dsn ; INFILE myfile; INPUT id $ &ptnames ; RUN ; %MEND readdata ; %readdata ; *---- DATA TEST commented out below is testefa.dta used ----- ; *---- used in Ferson and Rohlf's Fortran programs ----------- ; *---- This may be used in lieu of the data reading above ---- ; /* DATA TEST; INFILE CARDS; INPUT id $ x1 y1 x2 y2 x3 y3 x4 y4 x5 y5 x6 y6 x7 y7 x8 y8 x9 y9 x10 y10 x11 y11 x12 y12 ; CARDS ; 12 1 1 1 2 1 3 2 3 3 3 3 4 4 4 5 3 5 2 4 1 3 1 2 1 ; RUN ; *---------- Integer number of harmonics = harms ; */ PROC IML ; START xy2del(test,harms); USE &dsn ; READ ALL INTO r[COLNAME=coords ROWNAME=id]; k = (ncol(r))/2 ; /* PRINT k ; */ *--- Ensure that the number of harmonics is less than ---- ; *--- half the number of ordered pairs ---- ; harms = &harmonics ; IF &harmonics >= k/2 THEN DO ; harms = k/2; PRINT "harms reset to " ; PRINT harms ; END ; x = j(1,k); y = j(1,k); dx = j(1,k) ; dy = j(1,k) ; DO i=1 TO k ; x[,i]=r[,(i-1)*2+1]; y[,i]=r[,(i-1)*2+2]; END ; x1 = x[1,1]; y1 = y[1,1]; xmid=0.0 ; ymid=0.0 ; DO j=1 TO k-1 ; dx[1,j] = x[1,j+1] - x[1,j] ; dy[1,j] = y[1,j+1] - y[1,j] ; xmid=xmid+x[1,j]; ymid=ymid+y[1,j]; END ; dx[1,k] = x[1,1] - x[1,k] ; dy[1,k] = y[1,1] - y[1,k] ; xmid=(xmid+x[1,k])/k; ymid=(ymid+y[1,k])/k; *-- For location invariance ---- ; /* x1 = x1 - xmid ; y1 = y1 - ymid ; */ tsum = 0. ; xsum = 0. ; ysum = 0. ; a0=0. ; c0=0. ; xi = 0. ; yi = 0. ; t = j(1,k) ; rdx = j(1,k); rdy = j(1,k); DO j = 1 TO k ; t[1,j] = sqrt(dx[1,j]*dx[1,j] + dy[1,j]*dy[1,j]) ; rdx[1,j] = dx[1,j]/t[1,j] ; rdy[1,j] = dy[1,j]/t[1,j] ; tnew = tsum + t[1,j] ; IF (j > 1) THEN DO ; xi = xsum-rdx[1,j]*tsum ; yi = ysum-rdy[1,j]*tsum ; END ; t1 = t[1,j] ; t2 = tnew*tnew-tsum*tsum ; a0 = a0+.5*rdx[1,j]*t2+xi*t1 ; c0 = c0+.5*rdy[1,j]*t2+yi*t1 ; tsum = tnew ; xsum = xsum+dx[1,j] ; ysum = ysum+dy[1,j] ; END ; IF (tsum > 0.0) THEN DO ; a0 = x1+a0/tsum ; c0 = y1+c0/tsum ; END ; *---- LOOP ON HARMONICS (harms) ------- ; a = j(1,harms); b = j(1,harms) ; c = j(1,harms) ; d = j(1,harms) ; DO h=1 TO harms; fh = h ; fact1 = fh*6.283185/tsum ; angprv = 0.; asum = 0.; bsum = 0.; csum = 0.; dsum = 0.; *------ SUM OVER ALL POINTS ------------ ; DO j=1 TO k ; ang = angprv+t[1,j]*fact1 ; wtac = COS(ang)-COS(angprv); wtbd = SIN(ang)-SIN(angprv); angprv = ang ; asum = asum+rdx[1,j]*wtac ; bsum = bsum+rdx[1,j]*wtbd ; csum = csum+rdy[1,j]*wtac ; dsum = dsum+rdy[1,j]*wtbd; END ; fact2 = fact1*fh*3.141596 ; a[1,h] = asum/fact2 ; b[1,h] = bsum/fact2 ; c[1,h] = csum/fact2 ; d[1,h] = dsum/fact2 ; END ; *---------- ; invsize = 0 ; invrot = 0 ; invstart = 0 ; aa=a[1,1]; ba=b[1,1]; ca=c[1,1]; da=d[1,1]; denom=aa*aa+ca*ca-ba*ba-da*da ; theta=.5*ATAN2(2.*(aa*ba+ca*da),denom) ; r=cos(theta) ; s=sin(theta) ; astar=r*aa+s*ba ; cstar=r*ca+s*da ; psi=ATAN2(cstar,astar) ; estar=SQRT(astar*astar+cstar*cstar) ; IF (estar = 0.0) THEN estar=1.0 ; IF (invsize = 0) THEN estar=1.0 ; IF (invrot = 0) THEN psi=0.0 ; IF (invstart = 0) THEN theta=0.0 ; r=COS(psi); s=SIN(psi); a0n = (r*a0+s*c0)/estar ; c0n = (-s*a0+r*c0)/estar ; DO h = 1 TO harms ; thetan=theta*h ; rn=COS(thetan); sn=SIN(thetan); one=a[1,h]*r+c[1,h]*s ; two=b[1,h]*r+d[1,h]*s ; thr=c[1,h]*r-a[1,h]*s ; four=d[1,h]*r-b[1,h]*s ; a[1,h]=(rn*one+sn*two)/estar ; b[1,h]=(rn*two-sn*one)/estar ; c[1,h]=(rn*thr+sn*four)/estar ; d[1,h]=(rn*four-sn*thr)/estar ; end ; size=estar ; rotate=psi ; start=theta ; a0 = a0n ; c0 = c0n ; xmid=xmid+a0 ; ymid=ymid+c0 ; twopi = 6.283185 ; xe = j(1,k) ; ye = j(1,k) ; DO j = 1 TO k ; xi = a0 ; yi = c0 ; delang = twopi*(j-1.)/k ; ang = delang ; DO h = 1 TO harms ; ca = COS(ang); sa = SIN(ang) ; xi = xi+ca*a[1,h]+sa*b[1,h] ; yi = yi+ca*c[1,h]+sa*d[1,h] ; ang = ang+delang ; END ; Xe[1,j] = XI ; Ye[1,j] = YI ; END ; xname ={x y}; xtran = xe` ; ytran = ye`; newf = xtran||ytran ; PRINT newf ; CREATE mydata FROM newf [COLNAME=xname]; APPEND FROM newf ; print size ; print rotate ; print start ; print a0; print c0 ; print xmid ; print ymid ;print A0 ; print C0; print A[]; print b[]; print c[]; print d[]; print xe[]; print ye[]; FINISH xy2del ; RUN xy2del(test,&harmonics); QUIT ; PROC transpose DATA = &dsn OUT = transposed ; RUN ; DATA myx (RENAME=(col1=x) ) ; SET transposed ; IF _n_/2 -INT(_n_/2) NE 0 ; RUN ; DATA myy (RENAME=(col1=y)) ; SET transposed ; IF _n_/2 -INT(_n_/2) = 0 ; RUN ; DATA newf ; MERGE myx myy ; KEEP x y ; id = _N_ ; RUN ; QUIT ; DATA newf ; SET newf ; id = _N_ ; RUN ; DATA mydata (RENAME = (x = x1 y = y1)); SET mydata ; id = _N_; RUN ; DATA fin ; merge newf mydata ; BY id ; RUN ; DATA finlast ; SET fin ; IF _N_ = 1 ; RUN ; DATA fin ; SET fin finlast ; RUN ; TITLE1 &datasrc ; SYMBOL1 INTERPOL=JOIN COLOR = 'green' width=3 ; SYMBOL2 INTERPOL=JOIN COLOR = 'red' width=3; *TITLE2 "Original Data" ; TITLE2 "Dataset: &dsn Number of Points: &pairs " ; TITLE3 "Number of Harmonics: &harmonics " ; PROC GPLOT DATA = fin; PLOT y*x=1 y1*x1 = 2 /OVERLAY ; RUN ; QUIT ;