function coadpc2_mean(rmp)
%  a hacked version of coadpc2 to show movie morphing
%  along the vector connecting the population means
%
% COADPC2, Does Dynamic Principal Component Analysis for 2 Populations
%      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 coa2pm.m
%      Reads list of files from LISTFILESTR1 and LISTFILESTR2
% Inputs:
%     rmp - Robust Mean Parameter:
%             0 - (or not specified) use the ordinary sample mean
%             1 - Use the Spherical L1 location estimate
%             2 - Use the Elliptical (rescaled by MADs) L1 location estimate
% Outputs:
%     Only graphics, in current Figure

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


%  Set parameters and defaults according to number of input arguments
%
if nargin == 1 ;    %  1 argument input, use default
  irmp = 0 ;      %  default mean robustness parameter
else ;              %  was specified, use that
  irmp = rmp ;
end ;




%  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 IBASIS IRECENT COMAP NRPT NRIF ;
%  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
%  IBASIS - 1 for full (triangle), 2 for radial (3 row)
%  IRECENT - To recenter (or not) decompositions
%  COMAP - Color Map
%  NRPT - Number of times to Repeat Movie
%  NRIF - Number of times to Repeat Individual Movie Frames
          %    (controls speed of .mpeg version)

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
global FEATMPATH SUBDIRSTR LISTR ;
        %  path for storing Zernike feature matrices
        %  SUBDIRSTR - current subdirectory (of both indata and outputs)
        %  string for inputting list file
global LISTFILESTR1 LISTFILESTR2 SETLIST1PATH SETLIST2PATH ;
        %  list file strings
        %  strings to set paths to lists
global PC2STR PC2LOAD1 PC2LOAD2 ;
        %  PC2STR - String with title info for 2 population analyses
        %  PC2LOAD1 - 1st Population Loadings for 2 population analyses
        %  PC2LOAD2 - 2nd Population Loadings for 2 population analyses




%  Read in matrix of Zernike feature coefficients, Population 1
%
disp('      Reading Files for Population 1') ;
eval(SETLIST1PATH) ;
          %  sets LISTR to Full Path List, for Population 1
mfv1 = coaflr ;
d1 = size(mfv1,1) ;
n1 = size(mfv1,2) ;



%  Read in matrix of Zernike feature coefficients, Population 2
%
disp('      Reading Files for Population 2') ;
eval(SETLIST2PATH) ;
          %  sets LISTR to Full Path List, for Population 2
mfv2 = coaflr ;
d2 = size(mfv2,1) ;
n2 = size(mfv2,2) ;





%  Get mean feature vector
%
if irmp == 0 ;        %  then do original non-robust sample mean
  vmean1 = mean(mfv1')' ;
  vmean2 = mean(mfv2')' ;
  meanstr = 'mean' ;
elseif irmp == 1 ;    %  then do spherical robust sample mean
  vmean1 = rmean(mfv1')' ;
  vmean2 = rmean(mfv2')' ;
  meanstr = 'S. L1' ;
elseif irmp == 2 ;    %  then do elliptical robust sample mean
  vmad1 = gpmad(mfv1',0)' ;
          %  column vector of MADs of each row
  smfv1 = mfv1 .* vec2mat(1 ./ vmad1, n) ;
          %  make each coordinate have MAD 1
  vmean1 = rmean(smfv1')' ;
  vmean1 = vmean1 .* vmad1 ;

  vmad2 = gpmad(mfv2',0)' ;
          %  column vector of MADs of each row
  smfv2 = mfv2 .* vec2mat(1 ./ vmad2, n) ;
          %  make each coordinate have MAD 1
  vmean2 = rmean(smfv2')' ;
  vmean2 = vmean2 .* vmad2 ;
  meanstr = 'E. L1' ;
end ;



%  Combine into one big data set
if d1 == d2 ;
  d = d1 ;
else ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Error from coaspc2.m:          !!!') ;
  disp('!!!   Different Number of Zernikes   !!!') ;
  disp('!!!      for these 2 lists           !!!') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
end ;
mfv = [mfv1, mfv2] ;
n = n1 + n2 ;



%  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 


%  Create vector of frame orders
%
  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
vmorder = vec2mat(vmorder,NRIF) ;
vmorder = reshape(vmorder,1,size(vmorder,2)*NRIF) ;
          %  replicate frame number, NRIF times



%  Main Principal Component Block (combined data set)
%

%  Get mean feature vector
%
vmean = mean([vmean1'; vmean2'])' ;


%  Get residuals from mean
%
vsd = std(mfv')' ;
mresid = mfv - vec2mat(vmean,n) ;



dv = vmean1 - vmean ;
          %  vector from overall mean, towards mean
          %  of first population
dv = dv ./ sqrt(dv' * dv) ;
          %  make it a direction vector



vp = mresid' * dv ;
          %  vector of projections of data onto ev1
vprep = max(abs(vp - mean(vp))) ;
          %  "most representative pc", defined as MAD
mpv = dv * vprep ;
          %  "average vector" for display




%  Form mean +- representative eigenvectors
mp1p = vmean + mpv ;
mp1m = vmean - mpv ;





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

subplot(2,1,2) ;
  pccur = vp ;
          %  Matrix of current PC loadings
  pccur1 = pccur(1:n1) ;
          %  Current PC Loadings, for Population 1
  pccur2 = pccur(n1+1:n) ;
          %  Current PC Loadings, for Population 2
  [kde1, xgrid1] = gpkde(pccur1) ;
  [kde2, xgrid2] = gpkde(pccur2) ;
  plot(xgrid1,kde1,'r',xgrid2,kde2,'c') ;
    title(['Proj. on Mean Dif.']) ;
    vlim = axis ;
          %  Get axis values for drawing vertical lines
    left = -vprep ;
    right = vprep ;
    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([pccur2 pccur2],[bottom (top + bottom)/2], 'c') ;
      plot([pccur1 pccur1],[(top + bottom)/2 top], 'r') ;
          %  Show locations of projections 
          %  (pop'n 2 first, since likely more spread)
      plot(xgrid1,kde1,'r',xgrid2,kde2,'c') ;
          %  Do a replot to ensure seeing kde
      plot([vprep vprep], [bottom top], 'y') ;
          %  Show location of current projection
    hold off ;

mmov1 = [] ;
for iframe = 1:nframe ;
  w = (iframe - 1) / (nframe - 1) ;
  vfeat = (1 - w) * mp1p + w * mp1m ;
  linex = (1 - w) * vprep + w * (-vprep) ;
  mout = coazr(vfeat) ;

  subplot(2,1,1) ;
  pcolor(mxgrid,mygrid,mout) ;
    shading flat ;
    axis('square') ;
      caxis([0 1]) ;
    colormap(COMAP) ;
    hold on ;   % add centering type and bar showing time gone
      if IRECENT == 0 ;  %  then everything is centered at apex normals
        plot(0,0,'+k');
      elseif IRECENT == 1 ;  %  then everything is recentered
        plot(0,0,'ok');
      end ;

      barl = 0.4 * R0 ;
      bary = 0.95 * R0 ;
      barc = bary - barl / 2 ;
        if w < .5 ;
          plot([barc+(w-1/2)*barl, 1.05*barc],[bary,bary],'r') ;
        elseif w == .5 ;
          plot([barc, 1.05*barc],[bary,bary],'m') ;
        else ;
          plot([barc, 1.05*barc+(w-1/2)*barl],[bary,bary],'b') ;
        end ;
        vhand = get(gca,'Children') ;
        set(vhand(1),'LineWidth',5) ;
    hold off ;

    subplot(2,1,2) ;
      vchand = get(gca,'Children') ;

      set(vchand(1),'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 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) ;



%  Save 2 population analysis Loadings
%
mload = [pc1 pc2 pc3 pc4] ;
PC2LOAD1 = mload(1:n1,:) ;
PC2LOAD2 = mload(n1+1:n,:) ;
