function coadpca
% COADPCA, Does Dynamic Principal Component Analysis
%      Corneal Topography Analysis System
%   Creates a menu window for choosing movie, and creates
%   4 other windows, each with a kde and a movie
%   for doing Principal Component Analysis
%      Intended to be called by coapam.m
% Outputs:
%     Only graphics, in current Figure

%    Copyright (c) J. S. Marron 1998
%    Department of Statistics, University of North Carolina




%  Set local parameters
%
ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array
nframe = 21 ;
          %  number of movie frames
    


%  Set up global variables
%
%  First some general globals
global R0 NMAX COMAP NRPT ;
%  R0 - radius of analysis
          %    this gets reloaded when a feature vec is reloaded from file
%  NMAX - maximum number of Zernike columns
          %    this gets reloaded when a feature vec is reloaded from file
%  COMAP - Color Map
%  NRPT - Number of times to Repeat Movie
global SCREENWD SCREENHT MENUWD ;
        % overall sizes of display area, generic menu width
global VERTGAP HORIZGAP TITLHT BUTTONHT
        %  generic sizes of objects on menus
global DPCAHANDPOS DPCAM1POS DPCAM2POS DPCAM3POS DPCAM4POS ;
        %  window positioning vectors
global DPCAM1HAND DPCAM2HAND DPCAM3HAND DPCAM4HAND ;
        %  Handles for movie windows in Dynamic PC Analysis
global DPCAM1 DPCAM2 DPCAM3 DPCAM4 ;
        %  Movie matrices in Dynamic PC Analysis





%  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 

  imean = round((nframe + 1) / 2) ;
vmorder = imean:nframe ;
vmorder = [vmorder, ((nframe-1):-1:1)] ;
vmorder = [vmorder, (2:imean)] ;
          %  vector of frame orders to start and stop at mean



%  Read in matrix of Zernike feature coefficients
%
mfv = coaflr ;
d = size(mfv,1) ;
n = size(mfv,2) ;



%  Do Principal Components
%
sst = sum(sum(mfv.^2)) ;
%  get mean feature vector, and resids from that
vmean = mean(mfv')' ;
vsd = std(mfv')' ;
ssm = n * sum(vmean.^2) ;
ssmr = sst - ssm ;
mresid = mfv - vec2mat(vmean,n) ;

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

[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 ;

ev2 = meigvec(:,d-1) ;
          %  eigenvector with largest eigenvalue
pc2 = mresid1' * ev2 ;
          %  vector of projections of data onto ev2
pc2rep = max(abs(pc2 - mean(pc2))) ;
          %  "most representative pc", defined as MAD
mpc2 = ev2 * pc2rep ;
          %  "average vector" for display
mresid2 = (mresid1 - ev2 * pc2') ;
ss2r = sum(sum(mresid2.^2)) ;
ss2 = ss1r - ss2r ;

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 ;

ev4 = meigvec(:,d-3) ;
          %  eigenvector with largest eigenvalue
pc4 = mresid3' * ev4 ;
          %  vector of projections of data onto ev4
pc4rep = max(abs(pc4 - mean(pc4))) ;
          %  "most representative pc", defined as MAD
mpc4 = ev4 * pc4rep ;
          %  "average vector" for display
mresid4 = (mresid3 - ev4 * pc4') ;
ss4r = sum(sum(mresid4.^2)) ;
ss4 = ss3r - ss4r ;


%  Form mean +- representative eigenvectors
mpc1p = vmean + mpc1 ;
mpc1m = vmean - mpc1 ;
mpc2p = vmean + mpc2 ;
mpc2m = vmean - mpc2 ;
mpc3p = vmean + mpc3 ;
mpc3m = vmean - mpc3 ;
mpc4p = vmean + mpc4 ;
mpc4m = vmean - mpc4 ;





%  Create PC1 movie
%
DPCAM1HAND = figure ;
          %  create a new window
set(DPCAM1HAND, 'NumberTitle', 'off');
set(DPCAM1HAND,'Position', DPCAM1POS) ;
set(DPCAM1HAND,'Name','PC1') ;

subplot(2,1,2) ;
[kde, xgrid] = gpkde(pc1) ;
plot(xgrid,kde,'w') ;
  title('Proj 1') ;
  vlim = axis ;
          %  Get axis values for drawing vertical lines
  left = -pc1rep ;
  right = pc1rep ;
  range = right - left ;
  left = left - .05 * range ;
  right = right + .05 * range ;
  bottom = vlim(3) ;
  top = vlim(4) ;
  axis([left right bottom top]) ;
  hold on ;
    plot([pc1 pc1], [bottom (top+bottom)/2], 'c') ;
          %  Show locations of projections
    plot(xgrid,kde,'w') ;
          %  Do a replot to ensure seeing kde
    plot([pc1rep pc1rep], [bottom top], 'y') ;
          %  Show location of current projection
      tx = left + .1 * (right - left) ;
      ty = bottom + .9 * (top - bottom) ;
      ssstr = num2str(100 * ss1 / ssmr) ;
      tstr = [ssstr '% of MR'] ;
    text(tx,ty,tstr) ;
    title('Proj, PC1') ;
  hold off ;

mmov1 = [] ;
for iframe = 1:nframe ;
  w = (iframe - 1) / (nframe - 1) ;
  vfeat = (1 - w) * mpc1p + w * mpc1m ;
  linex = (1 - w) * pc1rep + w * (-pc1rep) ;
  mout = coazr(vfeat) ;

  subplot(2,1,1) ;
  pcolor(mxgrid,mygrid,mout) ;
    shading flat ;
    axis('square') ;
      caxis([0 1]) ;
    colormap(COMAP) ;
    hold on ;   % add bar showing time gone
      barl = 1.5 ;
      barc = 3.25 ;
      bary = 4 ;
        if w < .5 ;
          cstr = 'r' ;
        elseif w == .5 ;
          cstr = 'm' ;
        else ;
          cstr = 'b' ;
        end ;
      plot([barc,barc+(w-1/2)*barl],[bary,bary],cstr) ;
        vhand = get(gca,'Children') ;
        set(vhand(1),'LineWidth',5) ;
    hold off ;

    subplot(2,1,2) ;
      vchand = get(gca,'Children') ;
      set(vchand(2),'XData',[linex linex]) ;
          %  Move yellow line over to show location of current projection

   DPCAM1(:,iframe) = getframe(DPCAM1HAND) ;
end ;

DPCAM1 = DPCAM1(:,vmorder') ;
          %  reorder frames, to start and stop at mean

movie(DPCAM1HAND,DPCAM1,NRPT) ;








%  Create PC2 movie
%
DPCAM2HAND = figure ;
          %  create a new window
set(DPCAM2HAND, 'NumberTitle', 'off');
set(DPCAM2HAND,'Position', DPCAM2POS) ;
set(DPCAM2HAND,'Name','PC2') ;

subplot(2,1,2) ;
[kde, xgrid] = gpkde(pc2) ;
plot(xgrid,kde,'w') ;
  title('Proj 2') ;
  vlim = axis ;
          %  Get axis values for drawing vertical lines
  left = -pc2rep ;
  right = pc2rep ;
  range = right - left ;
  left = left - .05 * range ;
  right = right + .05 * range ;
  bottom = vlim(3) ;
  top = vlim(4) ;
  axis([left right bottom top]) ;
  hold on ;
    plot([pc2 pc2], [bottom (top+bottom)/2], 'c') ;
          %  Show locations of projections
    plot(xgrid,kde,'w') ;
          %  Do a replot to ensure seeing kde
    plot([pc1rep pc2rep], [bottom top], 'y') ;
          %  Show location of current projection
      tx = left + .1 * (right - left) ;
      ty = bottom + .9 * (top - bottom) ;
      ssstr = num2str(100 * ss2 / ssmr) ;
      tstr = [ssstr '% of MR'] ;
    text(tx,ty,tstr) ;
    title('Proj, PC2') ;
  hold off ;

mmov2 = [] ;
for iframe = 1:nframe ;
  w = (iframe - 1) / (nframe - 1) ;
  vfeat = (1 - w) * mpc2p + w * mpc2m ;
  linex = (1 - w) * pc2rep + w * (-pc2rep) ;
  mout = coazr(vfeat) ;

  subplot(2,1,1) ;
  colormap(COMAP) ;
  pcolor(mxgrid,mygrid,mout) ;
    shading flat ;
    axis('square') ;
      caxis([0 1]) ;
    colormap(COMAP) ;
    hold on ;   % add bar showing time gone
      barl = 1.5 ;
      barc = 3.25 ;
      bary = 4 ;
        if w < .5 ;
          cstr = 'r' ;
        elseif w == .5 ;
          cstr = 'm' ;
        else ;
          cstr = 'b' ;
        end ;
      plot([barc,barc+(w-1/2)*barl],[bary,bary],cstr) ;
        vhand = get(gca,'Children') ;
        set(vhand(1),'LineWidth',5) ;
    hold off ;

    subplot(2,1,2) ;
      vchand = get(gca,'Children') ;
      set(vchand(2),'XData',[linex linex]) ;
          %  Move yellow line over to show location of current projection

   DPCAM2(:,iframe) = getframe(DPCAM2HAND) ;
end ;

DPCAM2 = DPCAM2(:,vmorder') ;
          %  reorder frames, to start and stop at mean

movie(DPCAM2HAND,DPCAM2,NRPT) ;




%  Create PC3 movie
%
DPCAM3HAND = figure ;
          %  create a new window
set(DPCAM3HAND, 'NumberTitle', 'off');
set(DPCAM3HAND,'Position', DPCAM3POS) ;
set(DPCAM3HAND,'Name','PC3') ;

subplot(2,1,2) ;
[kde, xgrid] = gpkde(pc3) ;
plot(xgrid,kde,'w') ;
  title('Proj 3') ;
  vlim = axis ;
          %  Get axis values for drawing vertical lines
  left = -pc3rep ;
  right = pc3rep ;
  range = right - left ;
  left = left - .05 * range ;
  right = right + .05 * range ;
  bottom = vlim(3) ;
  top = vlim(4) ;
  axis([left right bottom top]) ;
  hold on ;
    plot([pc3 pc3], [bottom (top+bottom)/2], 'c') ;
          %  Show locations of projections
    plot(xgrid,kde,'w') ;
          %  Do a replot to ensure seeing kde
    plot([pc3rep pc3rep], [bottom top], 'y') ;
          %  Show location of current projection
      tx = left + .1 * (right - left) ;
      ty = bottom + .9 * (top - bottom) ;
      ssstr = num2str(100 * ss3 / ssmr) ;
      tstr = [ssstr '% of MR'] ;
    text(tx,ty,tstr) ;
    title('Proj, PC3') ;
  hold off ;

mmov3 = [] ;
for iframe = 1:nframe ;
  w = (iframe - 1) / (nframe - 1) ;
  vfeat = (1 - w) * mpc3p + w * mpc3m ;
  linex = (1 - w) * pc3rep + w * (-pc3rep) ;
  mout = coazr(vfeat) ;

  subplot(2,1,1) ;
  pcolor(mxgrid,mygrid,mout) ;
    shading flat ;
    axis('square') ;
      caxis([0 1]) ;
    colormap(COMAP) ;
    hold on ;   % add bar showing time gone
      barl = 1.5 ;
      barc = 3.25 ;
      bary = 4 ;
        if w < .5 ;
          cstr = 'r' ;
        elseif w == .5 ;
          cstr = 'm' ;
        else ;
          cstr = 'b' ;
        end ;
      plot([barc,barc+(w-1/2)*barl],[bary,bary],cstr) ;
        vhand = get(gca,'Children') ;
        set(vhand(1),'LineWidth',5) ;
    hold off ;

    subplot(2,1,2) ;
      vchand = get(gca,'Children') ;
      set(vchand(2),'XData',[linex linex]) ;
          %  Move yellow line over to show location of current projection

   DPCAM3(:,iframe) = getframe(DPCAM3HAND) ;
end ;

DPCAM3 = DPCAM3(:,vmorder') ;
          %  reorder frames, to start and stop at mean

movie(DPCAM3HAND,DPCAM3,NRPT) ;




%  Create PC4 movie
%
DPCAM4HAND = figure ;
          %  create a new window
set(DPCAM4HAND, 'NumberTitle', 'off');
set(DPCAM4HAND,'Position', DPCAM4POS) ;
set(DPCAM4HAND,'Name','PC4') ;

subplot(2,1,2) ;
[kde, xgrid] = gpkde(pc4) ;
plot(xgrid,kde,'w') ;
  title('Proj 4') ;
  vlim = axis ;
          %  Get axis values for drawing vertical lines
  left = -pc4rep ;
  right = pc4rep ;
  range = right - left ;
  left = left - .05 * range ;
  right = right + .05 * range ;
  bottom = vlim(3) ;
  top = vlim(4) ;
  axis([left right bottom top]) ;
  hold on ;
    plot([pc4 pc4], [bottom (top+bottom)/2], 'c') ;
          %  Show locations of projections
    plot(xgrid,kde,'w') ;
          %  Do a replot to ensure seeing kde
    plot([pc4rep pc4rep], [bottom top], 'y') ;
          %  Show location of current projection
      tx = left + .1 * (right - left) ;
      ty = bottom + .9 * (top - bottom) ;
      ssstr = num2str(100 * ss4 / ssmr) ;
      tstr = [ssstr '% of MR'] ;
    text(tx,ty,tstr) ;
    title('Proj, PC4') ;
  hold off ;

mmov4 = [] ;
for iframe = 1:nframe ;
  w = (iframe - 1) / (nframe - 1) ;
  vfeat = (1 - w) * mpc4p + w * mpc4m ;
  linex = (1 - w) * pc4rep + w * (-pc4rep) ;
  mout = coazr(vfeat) ;

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

    hold on ;   % add bar showing time gone
      barl = 1.5 ;
      barc = 3.25 ;
      bary = 4 ;
        if w < .5 ;
          cstr = 'r' ;
        elseif w == .5 ;
          cstr = 'm' ;
        else ;
          cstr = 'b' ;
        end ;
      plot([barc,barc+(w-1/2)*barl],[bary,bary],cstr) ;
        vhand = get(gca,'Children') ;
        set(vhand(1),'LineWidth',5) ;
    hold off ;

    subplot(2,1,2) ;
      vchand = get(gca,'Children') ;
      set(vchand(2),'XData',[linex linex]) ;
          %  Move yellow line over to show location of current projection

   DPCAM4(:,iframe) = getframe(DPCAM4HAND) ;
end ;

DPCAM4 = DPCAM4(:,vmorder') ;
          %  reorder frames, to start and stop at mean

movie(DPCAM4HAND,DPCAM4,NRPT) ;





%  Create Control window  (in new window)
%
dpcachand = figure ;
set(dpcachand, 'Position', DPCAHANDPOS);
set(dpcachand, 'NumberTitle', 'off');
set(dpcachand, 'Name', 'Push to Run Movie');



%  Create push button for PC 1 movie
%
buttony = VERTGAP ;
buttonwd = (MENUWD - 5 * HORIZGAP) / 4 ;
cbs1 = ['  figure(DPCAM1HAND) ;' ...
        '  colormap(COMAP) ;' ...
        '  movie(DPCAM1HAND,DPCAM1,NRPT) ;'] ;
          %  CallBack String 1: Runs movie for PC1
dpcam1hand = uicontrol(dpcachand,'Style','push',...
              'Position',[HORIZGAP buttony buttonwd BUTTONHT],...
              'String','PC1',...
              'CallBack',cbs1) ;



%  Create push button for PC 2 movie
%
buttony = VERTGAP ;
buttonwd = (MENUWD - 5 * HORIZGAP) / 4 ;
buttonx = buttonwd + 2 * HORIZGAP ;
cbs2 = ['  figure(DPCAM2HAND) ;' ...
        '  colormap(COMAP) ;' ...
        '  movie(DPCAM2HAND,DPCAM2,NRPT) ;'] ;
          %  CallBack String 2: Runs movie for PC2
dpcam2hand = uicontrol(dpcachand,'Style','push',...
              'Position',[buttonx buttony buttonwd BUTTONHT],...
              'String','PC2',...
              'CallBack',cbs2) ;



%  Create push button for PC 3 movie
%
buttony = VERTGAP ;
buttonwd = (MENUWD - 5 * HORIZGAP) / 4 ;
buttonx = 2 * buttonwd + 3 * HORIZGAP ;
cbs3 = ['  figure(DPCAM3HAND) ;' ...
        '  colormap(COMAP) ;' ...
        '  movie(DPCAM3HAND,DPCAM3,NRPT) ;'] ;
          %  CallBack String 3: Runs movie for PC3
dpcam3hand = uicontrol(dpcachand,'Style','push',...
              'Position',[buttonx buttony buttonwd BUTTONHT],...
              'String','PC3',...
              'CallBack',cbs3) ;



%  Create push button for PC 4 movie
%
buttony = VERTGAP ;
buttonwd = (MENUWD - 5 * HORIZGAP) / 4 ;
buttonx = 3 * buttonwd + 4 * HORIZGAP ;
cbs4 = ['  figure(DPCAM4HAND) ;' ...
        '  colormap(COMAP) ;' ...
        '  movie(DPCAM4HAND,DPCAM4,NRPT) ;'] ;
          %  CallBack String 4: Runs movie for PC4
dpcam4hand = uicontrol(dpcachand,'Style','push',...
              'Position',[buttonx buttony buttonwd BUTTONHT],...
              'String','PC4',...
              'CallBack',cbs4) ;

