disp('Running MATLAB script file lmtzcf35.m') ;
%
%    Makes Figure 3.5, 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 mean feature vector
%
vmean = mean(mfv')' ;


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

%  Get covariance matrix
%
mcov = cov(mresid') ;
          %  Get covariance matrix, transpose, since want "coordinates
          %                      as variables"

%  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
  
ev3 = meigvec(:,d-2) ;
          %  eigenvector with largest eigenvalue
pc3 = mresid2' * ev3 ;
          %  vector of projections of data onto ev3
pc3rep = max(abs(pc3 - mean(pc3))) ;
          %  "most representative pc", defined as MAD
mpc3 = ev3 * pc3rep ;
          %  "average vector" for display
mresid3 = (mresid2 - ev3 * pc3') ;
ss3r = sum(sum(mresid3.^2)) ;
ss3 = ss2r - ss3r ;


%  Form mean +- representative eigenvectors
mpc3p = vmean + mpc3 ;
mpc3m = vmean - mpc3 ;


%  Create images
%
imm = coazr(vmean) ;
impc3p = coazr(mpc3p) ;
impc3m = coazr(mpc3m) ;


%  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,impc3p)']) ;
  shading flat ;
  axis('square') ;
    caxis([0 1]) ;
  colormap(COMAP) ;

subplot(1,3,1) ;
eval(['pcolor(mxgrid,mygrid,impc3m)']) ;
  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\lmtzcf35.eps ;


