disp('Running MATLAB script file lmtzcf56.m') ;
%
%    Makes Figure 5.6, for the Locantore, Marron, Tripoli,
%    Zhang, Cohen paper.
%
%    lines mostly from coaspca.m



global R0 NMAX IBASIS IRECENT LISTR FEATMPATH SUBDIRSTR ;
          %  These are used by coaflr
global R0 NMAX IBASIS IRECENT INDATPATH HICURVE ;
          %  These are used by coazr

%  First define variables
R0 = 4.0 ;
NMAX = 10 ;
IBASIS = 1 ;
IRECENT = 0 ;
CRAD = 0 ;
CANG = 0 ;

FEATMPATH = '\matlab\cornean\zfm\' ;
SUBDIRSTR = 'milan\' ;
LISTFILESTR = 'prkpreo' ;
SETLISTPATH = 'LISTR = [FEATMPATH SUBDIRSTR LISTFILESTR ''.lst''] ; ' ;


%  Set up Keratron colormap
%
[COMAP, HICURVE] = coacmap ;


%  Set local parameters
%
ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array


%  Read in matrix of Zernike feature coefficients
%
eval(SETLISTPATH) ;
% sets LISTR (read-in file string)
mfv = coaflr ;
d = size(mfv,1) ;
n = size(mfv,2) ;



%  Calculate grid quantities
%
xgrid = linspace(-R0, R0, ngrid) ;
ygrid = xgrid ;	        %  Vectors of image coordinate values 
[mxgrid, mygrid] = meshgrid(xgrid,ygrid) ;
              		%  Matrices of image coordinate values 



%  Main Principal Component Block
%
sst = sum(sum(mfv.^2)) ;

%  Get elliptical mean feature vector
%
vmad = gpmad(mfv',0)' ;
          %  column vector of MADs of each row
smfv = mfv .* vec2mat(1 ./ vmad, n) ;
          %  make each coordinate have MAD 1
vmean = rmean(smfv')' ;
vmean = vmean .* vmad ;


%  Get residuals from mean
%
vsd = std(mfv')' ;
ssm = n * sum(vmean.^2) ;
ssmr = sst - ssm ;
mresid = mfv - vec2mat(vmean,n) ;

%  Get covariance matrix
%
vmad = gpmad(mresid',0)' ;
          %  column vector of MADs of each row
smresid = mresid .* vec2mat(1 ./ vmad, n) ;
          %  make each coordinate have MAD 1
vrad = sqrt(sum(smresid .^2)) ;
sphereresid = smresid' ./ vec2mat(vrad',d) ;
          %  Transpose, since want "coordinates as variables"
          %  vector of radii of each curve
          %  make each curve have "length" one
ellipsresid = sphereresid .* vec2mat(vmad', n) ;
          %  return each coordinate to original scale
mcov = cov(ellipsresid) ;
          %  Get covariance matrix of "ellipsed" data

%  Do eigen analysis
%
[meigvec,meigval] = eig(mcov) ;
veigval = (diag(meigval)) ;
          %  Get eigenvalues and eigenvectors

[veigval,vind] = sort(veigval) ;
meigvec = meigvec(:,vind) ;
          %  Sort to avoid possible ordering problems
  
ev1 = meigvec(:,d) ;
          %  eigenvector with largest eigenvalue
pc1 = mresid' * ev1 ;
          %  vector of projections of data onto ev1
pc1rep = max(abs(pc1 - mean(pc1))) ;
          %  "most representative pc", defined as MAD
mpc1 = ev1 * pc1rep ;
          %  "average vector" for display
mresid1 = (mresid - ev1 * pc1') ;
ss1r = sum(sum(mresid1.^2)) ;
ss1 = ssmr - ss1r ;


%  Form mean +- representative eigenvectors
mpc1p = vmean + mpc1 ;
mpc1m = vmean - mpc1 ;


%  Create images
%
imm = coazr(vmean) ;
impc1p = coazr(mpc1p) ;
impc1m = coazr(mpc1m) ;


%  Put down plots
%
%  Mean first
%
subplot(1,3,2) ;
pcolor(mxgrid,mygrid,imm) ;
  shading flat ;
  axis('square') ;
    caxis([0 1]) ;
  colormap(COMAP) ;


%  Now mean +- Pc's
%
subplot(1,3,3) ;
eval(['pcolor(mxgrid,mygrid,impc1p)']) ;
  shading flat ;
  axis('square') ;
    caxis([0 1]) ;
  colormap(COMAP) ;

subplot(1,3,1) ;
eval(['pcolor(mxgrid,mygrid,impc1m)']) ;
  shading flat ;
  axis('square') ;
    caxis([0 1]) ;
  colormap(COMAP) ;



  orient tall ;
  set(gcf,'PaperPosition',[1.5, 7.0, 5.5, 2.0]) ;
        %  1.5 inches on left and right, 5.5 interior, 8.5 total
        %  7.0 from bottom, 2.0 interior, 11 total
print -depsc \matlab\steve\ps\lmtzcf56.eps ;


