disp('Running MATLAB script file vire11.m') ;
%
%    Reads already stored Zernike Decomps, does PCA
%    for one person, and gives height differences between
%    mean and 1 (2) sd's out
%


nrgout = 26 ;
          %  number of radial gridpoints for output

nagout = 256 ;
          %  number of angular gridpoints for output





ipart = 1 ;      %  1 - PCA, and hieght difference maps
                 %  2 - Some Sum of Square Analyses

ipat = 0 ;       %  0 - loop through all
                 %  1 - Hickson
                 %  2 - Meadows
                 %  3 - Tesoro 
                 %  4 - Yi     


iR0 = 2 ;        %  1 - R0 = 4.2
                 %  2 - R0 = 3.0



npc = 20 ;       %  Number of Principal Components to Keep


irmp = 2 ;       %  0 - ordinary sample mean
                 %  1 - Spherical L1 Mean
                 %  2 - Elliptical L1 Mean

ircp = 2 ;       %  0 - ordinary covraince PCA
                 %  1 - Spherical PCA
                 %  2 - Elliptical PCA


%  Initialize parameters
%
if ipat == 0 ;     %  then loop through all patients
  ips = 1 ;
  ipe = 4 ;
else ;             %  then just do input patients
  ips = ipat ;
  ipe = ipat ;
end ;


global R0 NMAX IBASIS IRECENT INDATPATH HICURVE ;
NMAX = 10 ;
IBASIS = 1 ;
IRECENT = 0 ;


ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array
xgrid = linspace(-R0, R0, ngrid) ;
ygrid = xgrid ;	        %  Vectors of image coordinate values 


if iR0 == 1 ;
  outsuf = 'a' ;
  R0 = 4.2 ;
elseif iR0 == 2 ;
  outsuf = 'b' ;  
  R0 = 3.0 ;
end ;


INDATPATH = '\matlab\cornean\indata\' ;
          %  path for input Keratron files
[COMAP, HICURVE] = coacmap ;





if ipart == 1 ;    %  then do PCA, and get difference 
                   %                height maps


  %  Do appropriate file reads
  %
  for ip = ips:ipe ;

    if ip == 1 ;    
      patname = 'Hickson' ;
      featmpath = '\matlab\cornean\zfm\virehick\' ;
      listr = [featmpath 'hickall.lst'] ;
    elseif ip == 2 ;    
      patname = 'Meadows' ;
      featmpath = '\matlab\cornean\zfm\viremead\' ;
      listr = [featmpath 'meadall.lst'] ;
    elseif ip == 3 ;    
      patname = 'Tesoro' ;
      featmpath = '\matlab\cornean\zfm\vireteso\' ;
      listr = [featmpath 'tesoall.lst'] ;
    elseif ip == 4 ;    
      patname = 'Yi' ;
      featmpath = '\matlab\cornean\zfm\vireyi\' ;
      listr = [featmpath 'yiall.lst'] ;
    end ;



    %  read files, using lines from   coalwzr.m
    %
    fid = fopen(listr) ;
            %  assigns a file handle, (read only assumed)

    if (fid == -1) ;      %  then file not found, so report
      disp('   ') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      disp('!!!   Warning from vire2.m:   can''t find file:') ;
      disp(['!!!       ' listr]) ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    end ;


    %  Read in list of filenames
    %
    minfilestr = fscanf(fid,'%s',[8,inf]) ;
              %  formatted read, assumes single column
    fclose(fid) ;

    minfilestr = minfilestr' ;
            %  each row corresponds to one filename


    %  loop through files (one for each image)
    %
    mbetahat = [] ;
    ie = size(minfilestr,1) ;
    for i = 1:ie ;

      infilestr = minfilestr(i,:) ;

      disp(['  Loading ' num2str(i) ' of ' num2str(ie) ...
                               ' for ' patname ', file: ' infilestr ]) ;


      zmistr = [featmpath infilestr '.mat'] ;
            %  full path and filename 
      eval(['load ' zmistr]) ;
            %  loads this Zernike feature matrix and parameter values

      mbetahat = [mbetahat, betahat] ;



    end ;







    %  Do PCA
    %
    d = size(mbetahat,1) ;
    n = size(mbetahat,2) ;


    %  Get mean feature vector
    %
    if irmp == 0 ;        %  then do original non-robust sample mean
      vmean = mean(mbetahat')' ;
      meanstr = 'mean' ;
    elseif irmp == 1 ;    %  then do spherical robust sample mean
      vmean = rmean(mbetahat')' ;
      meanstr = 'S. L1' ;
    elseif irmp == 2 ;    %  then do elliptical robust sample mean
      vmad = gpmad(mbetahat',0)' ;
            %  column vector of MADs of each row
      smbetahat = mbetahat .* vec2mat(1 ./ vmad, n) ;
            %  make each coordinate have MAD 1
      vmean = rmean(smbetahat')' ;
      vmean = vmean .* vmad ;
      meanstr = 'E. L1' ;
    end ;

    %  Get residuals from mean
    %
    mresid = mbetahat - vec2mat(vmean,n) ;

    %  Get covariance matrix
    %
    if ircp == 0 ;        %  then do original non-robust covariance matrix
      mcov = cov(mresid') ;
          %  Get covariance matrix, transpose, since want "coordinates
          %                      as variables"
      covstr = 'PC' ;
    elseif ircp == 1 ;    %  then do spherical robust sample mean
      vrad = sqrt(sum(mresid .^2)) ;
          %  Transpose, since want "coordinates as variables"
          %  vector of radii of each curve
      sphereresid = mresid' ./ vec2mat(vrad',d) ;
          %  make each curve have "length" one
      mcov = cov(sphereresid) ;
          %  Get covariance matrix of "sphered" data
      covstr = 'SPC' ;
    elseif ircp == 2 ;    %  then do elliptical robust sample mean
      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
      covstr = 'EPC' ;
    end ;


    %  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
  

    mev = meigvec(:,(d-1):d) ;
          %  matrix of eigenvectors with 2 largest eigenvalues
    mpc = mresid' * mev ;
          %  matrix of projections of data onto 
          %  eigenvectors with largest eigenvalues
    mpc = mpc(:,2:-1:1) ;
          %  reverse the columns, to put the projections
          %  corresponding to the largest eigenvalues on the left


    pcsd = std(mpc) ;
          %  standard deviations of projections


    mm1sdpc1 = pcsd(1) * meigvec(:,1) ;
          %  Zernike coefficents, for mean minus 1 sd, for pc1


    %  Reconstruct, using lines modified from coazr.m:
%    %  Construct matrix of curvatures
%%
%if IBASIS == 1 ;    %  then use full basis
%  [mh,mhp,mhpp] = zrecon(featv,NMAX,mrgrid/R0,mthetag,R0) ;
%			%  creates vectors of fitted values
%			%  for a given data set (ht, deriv, 2deriv)
%elseif IBASIS == 2 ;   %  then use radial basis
%  [mh,mhp,mhpp] = rrecon(featv,NMAX,mrgrid/R0,mthetag,R0) ;
%			%  creates vectors of fitted values
%			%  for a given data set (ht, deriv, 2deriv)
%end ;



    vrgrid = linspace(0,R0,nrgout) ;
    vagrid = linspace(0,2*pi,nagout)' ;

    mrgrid = vec2mat(vrgrid,nagout) ;
    magrid = vec2mat(vagrid,nrgout) ;

    mrgrid = reshape(mrgrid,nrgout*nagout,1) ;
    magrid = reshape(magrid,nrgout*nagout,1) ;
          %  turn into big column vectors

    [mh,mhp,mhpp] = zrecon(mm1sdpc1,NMAX,mrgrid/R0,magrid,R0) ;


    mrgrid = reshape(mrgrid,nagout,nrgout) ;
    magrid = reshape(magrid,nagout,nrgout) ;
    mh = reshape(mh,nagout,nrgout) ;




    surf(mrgrid .* cos(magrid),mrgrid .* sin(magrid),mh) ;
      title(patname) ;
      shading flat ;
      colormap(jet) ;   
      view(0,90) ;
      colorbar ;

    gppause ;

    surfl(mrgrid .* cos(magrid),mrgrid .* sin(magrid),mh) ;
      title(patname) ;
      shading flat ;
      view(45,45) ;
      colormap(gray) ;   


    gppause ;


  end ;


end ;


