disp('MATLAB Script File:  cosd21.m') ;
%
%   For the Cornea Study, this does a "surfaces as data" approach
%       3rd Generation, using complete data sets, 
%           better eignevector rpresentation, and movie versions
%   This does standard principle components 
%
format compact ;


igroup = 20    %  1 - Repeatability Study (all)
              %  20 - PRK Study (Pre Op)
              %  21 - PRK Study (1 mo)
              %  26 - PRK Study (6 mo)
              %  212 - PRK Study (12 mo)
              %  300 - Kerataconic (normal)
              %  301 - Kerataconic (susp map)
              %  302 - Kerataconic (unilat)
              %  303 - Kerataconic (clinical)

imake = 0 ;      %  0 - do everything
                 %  1 - do eigebanalysis, make figure 1, mean and kde plots
                 %  2 - make figures with movies
                 %  3 - make GUI menu, to run movies in fig. 3
                 %  4 - make mpeg files

              
isurf = 3 ;   %  1 - for height map
              %  2 - for height map - sphere
              %  3 - for curvature map

ibas = 1 ;    %  1 - Radial Fourier - Legendre
              %  2 - Diameter Fourier - Legendre
              %            (not working yet)
              %  3 - Zernike
              %            (not working yet)



%  Set up Figure windows
%
fhand1 = figure(1) ;
  set(fhand1,'Position',[10,370,620,130]) ;

fhand2 = figure(2) ;
  set(fhand2,'Position',[120,100,130,130]) ;

fhand3 = figure(3) ;
  set(fhand3,'Position',[250,100,130,130]) ;

fhand4 = figure(4) ;
  set(fhand4,'Position',[380,100,130,130]) ;

fhand5 = figure(5) ;
  set(fhand5,'Position',[510,100,130,130]) ;

fhand6 = figure(6) ;
  set(fhand6,'Position',[5,100,110,130]) ;



%  Set file save string
%
fsavstr = ['cs21p' num2str(igroup)] ;





if imake == 1  | imake == 0 ;

  r0 = 4.2 ;    %  Usual radius of analysis

  ngrid = 64 ; %  Construct surfaces as ngrid x ngrid array

  xgrid = linspace(-r0,r0,ngrid) ;
  ygrid = xgrid ;
              %  Vectors of image coordinate values 

  [mxgrid, mygrid] = meshgrid(xgrid,ygrid) ;
              %  Matrices of image coordinate values 
              
  mrgrid = sqrt(mxgrid.^2 + mygrid.^2) ;
              %  Matrices of image polar coordinate radii
  mthetag = atan2(mxgrid,mygrid) ;
              %  Matrices of image polar coordinate angles (-pi,pi]
    vind = find(mthetag < 0) ;
  mthetag(vind) = mthetag(vind) + 2 * pi ;
              %  Matrices of image polar coordinate angles [0,2pi)
          
  mcrop = ones(ngrid,ngrid) ;            
  vind = find(mrgrid > r0) ;
              %  Vector of indices where image is outside circle
  mcrop(vind) = nan * ones(length(vind),1) ;
              %  Cropping Matrix, multiply by this to cut off images
              %  outside circle of radius r_0.       

  nr = 10 ;
  nc = 7 ;


  if igroup == 1 ;  %  Work with Repeatability population
    imagstr = 'Rep''t Pop''n' ;
    disp(['Loading feature matrices for ' imagstr]) ;
    %  First load in file info from file rptkey.dat, 
    %  which was modified from Jinting's MATLAB script normind.m
    %  (everything but the data was edited out).
    %    1st col is patient id
    %    2nd col is file id (add last 4 digits to "rpfm" to get file suffix)
    %    3rd col is "eye"
    %    4th col is pupil distance
    %    5th col is pupil direction
    fid = fopen('\matlab\steve\cornea\fm\rptkey.dat','r') ;
      frmtstr = '%g %g %g %g %g' ;
        mfilein = fscanf(fid,frmtstr,[5,inf]) ;
        mfilein = mfilein' ;
            %  since data rows are read in as columns
      fclose(fid) ;

    nimag = size(mfilein,1) ;
    vpatind = mfilein(:,1) ;
            %  vector of patient id's
    vfileid = mfilein(:,2) ;
            %  vector of file id's
    veyeind = mfilein(:,3) ;
            %  vector of eye indices
    mfilein = 0 ;
            %  won't be using this again


    %  Loop through feature matrix files, and combine into big data matrix
    mfv = [] ;
    for iimag = 1:nimag ;
      eval(['load \matlab\steve\cornea\fm\rpfm' ...
                               num2str(vfileid(iimag)) '.mat']) ;
            %  This load Jinting's previously saved stuff
            %  which is a number of variables, especially fm,
            %  a 19 x 7 "feature matrix"
      mfv = [mfv fm(:)] ;
            %  put matrix into a column vector and add to right
            %  Including the feature matrix for this image
    end ;


  elseif igroup < 300 ;   %  work with prk population
    disp('Working on PRK Population ') ;
    %  First load in file info from file prkkey.dat, 
    %  which was modified from fcfilesn.txt
    %  (everything but the data was edited out).
    %    1st col is patient id
    %    2nd col is Time (0,1,6,12 months)
    %    3nd col is file id (add last 4 digits to "fm42" to get file suffix)
    %    4th col is "eye"
    %    5th col is pupil distance
    %    6th col is pupil direction
    fid = fopen('\matlab\steve\cornea\fm\prkkey.dat','r') ;
      frmtstr = '%g %g %g %g %g %g' ;
        mfilein = fscanf(fid,frmtstr,[6,inf]) ;
        mfilein = mfilein' ;
            %  since data rows are read in as columns
      fclose(fid) ;

    vpatind = mfilein(:,1) ;
            %  vector of patient id's
    vtime = mfilein(:,2) ;
            %  vector of time's
    vfileid = mfilein(:,3) ;
            %  vector of file id's
    veyeind = mfilein(:,4) ;
            %  vector of eye indices
    mfilein = 0 ;
            %  won't be using this again


    if igroup == 20 ;
      imagstr = 'PRK Pre-op' ;
      disp(['Loading feature matrices for ' imagstr]) ;
      vflag = (vtime == 0) ;
    elseif igroup == 21 ;
      imagstr = 'PRK 1 Mo' ;
      disp(['Loading feature matrices for ' imagstr]) ;
      vflag = (vtime == 1) ;
    elseif igroup == 26 ;
      imagstr = 'PRK 6 Mo' ;
      disp(['Loading feature matrices for ' imagstr]) ;
      vflag = (vtime == 6) ;
    elseif igroup == 212 ;
      imagstr = 'PRK 12 Mo' ;
      disp(['Loading feature matrices for ' imagstr]) ;
      vflag = (vtime == 12) ;
    end ;

    vpatind = vpatind(vflag) ;
    vfileid = vfileid(vflag) ;
    veyeind = veyeind(vflag) ;
    nimag = length(vpatind) ;


    %  Loop through feature matrix files, and combine into big data matrix
    mfv = [] ;
    for iimag = 1:nimag ;
      fid = vfileid(iimag) ;
      if fid < 100 ;
        fids = ['00' num2str(fid)] ;
      elseif fid < 1000 ;
        fids = ['0' num2str(fid)] ;
      else ;
        fids = num2str(fid) ;
      end ;
      eval(['load \matlab\steve\cornea\fm\fm42' fids '.mat']) ;
            %  This load Jinting's previously saved stuff
            %  which is a number of variables, especially fm,
            %  a 19 x 7 "feature matrix"
      mfv = [mfv fm(:)] ;
            %  put matrix into a column vector and add to right
            %  Including the feature matrix for this image
    end ;


  else ;    %  work with Keratoconus group
    disp('Working on Ker Population ') ;
    %  First load in file info from file kerkey.dat, 
    %  which was modified from fcfilesn.txt
    %  (everything but the data was edited out).
    %    1st col is patient id
    %    2nd col is Ker type Code
    %    3nd col is file id (add last 4 digits to "fm42" to get file suffix)
    %    4th col is "eye"
    %    5th col is pupil distance
    %    6th col is pupil direction
    fid = fopen('\matlab\steve\cornea\fm\kerkey.dat','r') ;
      frmtstr = '%g %g %g %g %g %g' ;
        mfilein = fscanf(fid,frmtstr,[6,inf]) ;
        mfilein = mfilein' ;
              %  since data rows are read in as columns
      fclose(fid) ;

    vpatind = mfilein(:,1) ;
            %  vector of patient id's
    vkercode = mfilein(:,2) ;
            %  vector of time's
    vfileid = mfilein(:,3) ;
            %  vector of file id's
    veyeind = mfilein(:,4) ;
            %  vector of eye indices
    mfilein = 0 ;
            %  won't be using this again



    if igroup == 300 ;
      imagstr = 'Ker Normal' ;
      disp(['Loading feature matrices for ' imagstr]) ;
      vflag = (vkercode == 0) ;
    elseif igroup == 301 ;
      imagstr = 'Ker sus map' ;
      disp(['Loading feature matrices for ' imagstr]) ;
      vflag = (vkercode == 1) ;
    elseif igroup == 302 ;
      imagstr = 'Ker unilat' ;
      disp(['Loading feature matrices for ' imagstr]) ;
      vflag = (vkercode == 2) ;
    elseif igroup == 303 ;
      imagstr = 'Ker clinical' ;
      disp(['Loading feature matrices for ' imagstr]) ;
      vflag = (vkercode == 3) ;
    end ;

    vpatind = vpatind(vflag) ;
    vfileid = vfileid(vflag) ;
    veyeind = veyeind(vflag) ;
    nimag = length(vpatind) ;


    %  Loop through feature matrix files, and combine into big data matrix
    mfv = [] ;
    for iimag = 1:nimag ;
      fid = vfileid(iimag) ;
      if fid < 100 ;
        fids = ['00' num2str(fid)] ;
      elseif fid < 1000 ;
        fids = ['0' num2str(fid)] ;
      else ;
        fids = num2str(fid) ;
      end ;
      eval(['load \matlab\steve\cornea\fm\fm42' fids '.mat']) ;
              %  This load Jinting's previously saved stuff
              %  which is a number of variables, especially fm,
              %  a 19 x 7 "feature matrix"
      mfv = [mfv fm(:)] ;
              %  put matrix into a column vector and add to right
              %  Including the feature matrix for this image
    end ;


  end ;

  fclose('all') ;




  disp(['feature matrix has size ' num2str(size(mfv,1)) ' rows ']) ;
  disp(['and ' num2str(size(mfv,2)) ' cols']) ;





  disp('Working on eigenector decomposition') ;
  d = size(mfv,1) ;
  n = size(mfv,2) ;
  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 = median(abs(pc1 - mean(pc1))) ;
  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 = median(abs(pc2 - mean(pc2))) ;
  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 = median(abs(pc3 - mean(pc3))) ;
  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 = median(abs(pc4 - mean(pc4))) ;
  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 ;


  %  Reshape representative guys to feature matrix form
  fmmean = reshape(vmean,19,7) ;
  fmpc1p = reshape(mpc1p,19,7) ;
  fmpc2p = reshape(mpc2p,19,7) ;
  fmpc3p = reshape(mpc3p,19,7) ;
  fmpc4p = reshape(mpc4p,19,7) ;
  fmpc1m = reshape(mpc1m,19,7) ;
  fmpc2m = reshape(mpc2m,19,7) ;
  fmpc3m = reshape(mpc3m,19,7) ;
  fmpc4m = reshape(mpc4m,19,7) ;
  

  %  start plotting stuff, from corecon5.m
  disp('Generating surface') ; 
  if isurf == 1 | isurf == 2 ;   %  then work with height map
    surfstr = ' H''ts' ;

    if ibas == 1 ;     %  then use the radial Fourier - Legendre basis

      moutm = corh(fmmean,mrgrid,mthetag,r0) ;
      mout1p = corh(fmpc1p,mrgrid,mthetag,r0) ;
      mout2p = corh(fmpc2p,mrgrid,mthetag,r0) ;
      mout3p = corh(fmpc3p,mrgrid,mthetag,r0) ;
      mout4p = corh(fmpc4p,mrgrid,mthetag,r0) ;
      mout1m = corh(fmpc1m,mrgrid,mthetag,r0) ;
      mout2m = corh(fmpc2m,mrgrid,mthetag,r0) ;
      mout3m = corh(fmpc3m,mrgrid,mthetag,r0) ;
      mout4m = corh(fmpc4m,mrgrid,mthetag,r0) ;
          %  Radial Fourier - Legendre reconstructions of Heights

    elseif ibas == 2 ;     %  then use the diameter Fourier - Legendre basis

    elseif ibas == 3 ;     %  then use the Zernike basis

    end ;

    if isurf == 2 ;   %  Then need to subtract sphere
      surfstr = [surfstr ' - Sphere'] ;

  %    minmout = min(min(mout))
  %    maxmout = max(max(mout))
            %  use these lines to decide which sphere to subtract

      sphht = spherad - sqrt(spherad.^2 - mrgrid.^2) ;

  %    minsph = min(min(sphht))
  %    maxsph = max(max(sphht))
  %    pause
            %  use these lines to decide which sphere to subtract


      moutm = moutm + sphht ;
      mout1p = mout1p + sphht ;
      mout2p = mout2p + sphht ;
      mout3p = mout3p + sphht ;
      mout4p = mout4p + sphht ;
      mout1m = mout1m + sphht ;
      mout2m = mout2m + sphht ;
      mout3m = mout3m + sphht ;
      mout4m = mout4m + sphht ;
            %  careful about whether this is + or -
    end ;

  elseif isurf == 3 ;   %  then work with curvature map
    surfstr = ' Curv''s' ;

    if ibas == 1 ;     %  then use the radial Fourier - Legendre basis

      [mhm,mhpm,mhppm] = corhpp(fmmean',mrgrid,mthetag,r0) ;
      [mh1p,mhp1p,mhpp1p] = corhpp(fmpc1p',mrgrid,mthetag,r0) ;
      [mh2p,mhp2p,mhpp2p] = corhpp(fmpc2p',mrgrid,mthetag,r0) ;
      [mh3p,mhp3p,mhpp3p] = corhpp(fmpc3p',mrgrid,mthetag,r0) ;
      [mh4p,mhp4p,mhpp4p] = corhpp(fmpc4p',mrgrid,mthetag,r0) ;
      [mh1m,mhp1m,mhpp1m] = corhpp(fmpc1m',mrgrid,mthetag,r0) ;
      [mh2m,mhp2m,mhpp2m] = corhpp(fmpc2m',mrgrid,mthetag,r0) ;
      [mh3m,mhp3m,mhpp3m] = corhpp(fmpc3m',mrgrid,mthetag,r0) ;
      [mh4m,mhp4m,mhpp4m] = corhpp(fmpc4m',mrgrid,mthetag,r0) ;

    elseif ibas == 2 ;     %  then use the diameter Fourier - Legendre basis

    elseif ibas == 3 ;     %  then use the Zernike basis
    end ;

    moutm = -mhppm ./ (1 + mhpm.^2).^(3/2) ;
    mout1p = -mhpp1p ./ (1 + mhp1p.^2).^(3/2) ;
    mout2p = -mhpp2p ./ (1 + mhp2p.^2).^(3/2) ;
    mout3p = -mhpp3p ./ (1 + mhp3p.^2).^(3/2) ;
    mout4p = -mhpp4p ./ (1 + mhp4p.^2).^(3/2) ;
    mout1m = -mhpp1m ./ (1 + mhp1m.^2).^(3/2) ;
    mout2m = -mhpp2m ./ (1 + mhp2m.^2).^(3/2) ;
    mout3m = -mhpp3m ./ (1 + mhp3m.^2).^(3/2) ;
    mout4m = -mhpp4m ./ (1 + mhp4m.^2).^(3/2) ;
    moutm = moutm * 337.5 ;
    mout1p = mout1p * 337.5 ;
    mout2p = mout2p * 337.5 ;
    mout3p = mout3p * 337.5 ;
    mout4p = mout4p * 337.5 ;
    mout1m = mout1m * 337.5 ;
    mout2m = mout2m * 337.5 ;
    mout3m = mout3m * 337.5 ;
    mout4m = mout4m * 337.5 ;
          %  put on scale of diopters

  end ;




  %  do explicit flip of x's and y's, since pcolor wants:
  %            rows as y's,   cols as x's
  moutm = moutm' ;
  mout1p = mout1p' ;
  mout2p = mout2p' ;
  mout3p = mout3p' ;
  mout4p = mout4p' ;
  mout1m = mout1m' ;
  mout2m = mout2m' ;
  mout3m = mout3m' ;
  mout4m = mout4m' ;







  if isurf == 2 ;
          %  Then setup Keratron height color map

    %  These use the "full range" Keratron data, from the email message:
    %             cocolora.msg
    %   which were copied into the ASCII file 
    %             cocolor1.dat
    %  for importing into matlab programs
    %    locurve & hicurve are upper and lower ends of the "range"
    %    for each color
    load cornea\cocolor1.dat ;
    cyan = cocolor1(:,3) ;
    magenta = cocolor1(:,4) ;
    yellow = cocolor1(:,5) ;
    black = cocolor1(:,6) ;
    %  Don't use columns 1 and 2, because those are for curvature

    red   = 1 - cyan + black ;
    green = 1 - magenta + black;
    blue  = 1 - yellow + black ;

    comap = [red green blue] ;
    ncomap = size(comap,1) ;
          %  number of rows in the color map ;

    %  make sure values are positive
    %     (since twiddling "black" can leave this range)
    flag = comap < 0 ;
    fflag = find(flag) ;
    comap(fflag) = zeros(length(fflag),1) ;

    %  make sure values are less than 1
    %     (since twiddling "black" can leave this range)
    flag = comap > 1 ;
    fflag = find(flag) ;
    comap(fflag) = ones(length(fflag),1) ;

    colormap(comap) ;

        cstep = cstep / 1000 ;
          %  cstep is "step size in color scale", different for
          %  each image

        u01 = linspace(0,1,ncomap+1) ;
          lbd = -11.5 * cstep ;
          ubd = lbd + ncomap * cstep ;
        colgrid = linspace(lbd,ubd,ncomap+1) ;
          %  equally spaced grid on desired color scale

      vmout = mout(:) ;
          %  vector version for easier manipulation
          %  NOTE: interp1 seems to need vector argument
        indlo = find(vmout < lbd) ;
      vmout(indlo) = lbd * ones(length(indlo),1) ;
        indup = find(vmout > ubd) ;
      vmout(indup) = lbd * ones(length(indup),1) ;
          %  NOTE: interp1 seems to need values within grid

      colcurv = interp1(colgrid, u01, vmout) ;
          %  use this to determine vector of colors representing
          %  curvature, when plotting.  This transforms curvature 
          %  to get Keratron colors, when using comap as colormap.
    moutm = reshape(colcurv,size(moutm)) ;  
    mout1p = reshape(colcurv,size(mout1p)) ;  
    mout2p = reshape(colcurv,size(mout2p)) ;  
    mout3p = reshape(colcurv,size(mout3p)) ;  
    mout4p = reshape(colcurv,size(mout4p)) ;  
    mout1m = reshape(colcurv,size(mout1m)) ;  
    mout2m = reshape(colcurv,size(mout2m)) ;  
    mout3m = reshape(colcurv,size(mout3m)) ;  
    mout4m = reshape(colcurv,size(mout4m)) ;  

  elseif isurf == 3 ;
          %  Then setup Keratron color map

    %  These use the "full range" Keratron data, from the email message:
    %             cocolora.msg
    %   which were copied into the ASCII file 
    %             cocolor1.dat
    %  for importing into matlab programs
    %    locurve & hicurve are upper and lower ends of the "range"
    %    for each color
    load \matlab\steve\cornea\cocolor1.dat ;
    locurve = cocolor1(:,1) ;
    hicurve = cocolor1(:,2) ;
    cyan = cocolor1(:,3) ;
    magenta = cocolor1(:,4) ;
    yellow = cocolor1(:,5) ;
    black = cocolor1(:,6) ;

    red   = 1 - cyan + black ;
    green = 1 - magenta + black;
    blue  = 1 - yellow + black ;

    comap = [red green blue] ;
    ncomap = size(comap,1) ;
          %  number of rows in the color map ;

    %  make sure values are positive
    flag = comap < 0 ;
    fflag = find(flag) ;
    comap(fflag) = zeros(length(fflag),1) ;

    %  make sure values are less than 1
    flag = comap > 1 ;
    fflag = find(flag) ;
    comap(fflag) = ones(length(fflag),1) ;

      vmoutm = moutm(:) ;
      vmout1p = mout1p(:) ;
      vmout2p = mout2p(:) ;
      vmout3p = mout3p(:) ;
      vmout4p = mout4p(:) ;
      vmout1m = mout1m(:) ;
      vmout2m = mout2m(:) ;
      vmout3m = mout3m(:) ;
      vmout4m = mout4m(:) ;
          %  vector version for easier manipulation
          %  NOTE: interp1 seems to need vector argument
        indlo = find(vmoutm < 0) ;
      vmoutm(indlo) = 0 * ones(length(indlo),1) ;
        indup = find(vmoutm > max(hicurve)) ;
      vmoutm(indup) = max(hicurve) * ones(length(indup),1) ;
          %  NOTE: interp1 seems to need values within grid
        indlo = find(vmout1p < 0) ;
      vmout1p(indlo) = 0 * ones(length(indlo),1) ;
        indup = find(vmout1p > max(hicurve)) ;
      vmout1p(indup) = max(hicurve) * ones(length(indup),1) ;
        indlo = find(vmout2p < 0) ;
      vmout2p(indlo) = 0 * ones(length(indlo),1) ;
        indup = find(vmout2p > max(hicurve)) ;
      vmout2p(indup) = max(hicurve) * ones(length(indup),1) ;
        indlo = find(vmout3p < 0) ;
      vmout3p(indlo) = 0 * ones(length(indlo),1) ;
        indup = find(vmout3p > max(hicurve)) ;
      vmout3p(indup) = max(hicurve) * ones(length(indup),1) ;
        indlo = find(vmout4p < 0) ;
      vmout4p(indlo) = 0 * ones(length(indlo),1) ;
        indup = find(vmout4p > max(hicurve)) ;
      vmout4p(indup) = max(hicurve) * ones(length(indup),1) ;
        indlo = find(vmout1m < 0) ;
      vmout1m(indlo) = 0 * ones(length(indlo),1) ;
        indup = find(vmout1m > max(hicurve)) ;
      vmout1m(indup) = max(hicurve) * ones(length(indup),1) ;
        indlo = find(vmout2m < 0) ;
      vmout2m(indlo) = 0 * ones(length(indlo),1) ;
        indup = find(vmout2m > max(hicurve)) ;
      vmout2m(indup) = max(hicurve) * ones(length(indup),1) ;
        indlo = find(vmout3m < 0) ;
      vmout3m(indlo) = 0 * ones(length(indlo),1) ;
        indup = find(vmout3m > max(hicurve)) ;
      vmout3m(indup) = max(hicurve) * ones(length(indup),1) ;
        indlo = find(vmout4m < 0) ;
      vmout4m(indlo) = 0 * ones(length(indlo),1) ;
        indup = find(vmout4m > max(hicurve)) ;
      vmout4m(indup) = max(hicurve) * ones(length(indup),1) ;

    colormap(comap) ;

      u01 = linspace(0,1,ncomap+1) ;
      colcurv = interp1([0; hicurve], u01, vmoutm) ;
          %  use this to determine vector of colors representing
          %  curvature, when plotting.  This transforms curvature 
          %  to get Keratron colors, when using comap as colormap.
          %  NOTE: interp1 seems to need vector argument
    moutm = reshape(colcurv,size(moutm)) ;  
      colcurv = interp1([0; hicurve], u01, vmout1p) ;
    mout1p = reshape(colcurv,size(mout1p)) ;  
      colcurv = interp1([0; hicurve], u01, vmout2p) ;
    mout2p = reshape(colcurv,size(mout2p)) ;  
      colcurv = interp1([0; hicurve], u01, vmout3p) ;
    mout3p = reshape(colcurv,size(mout3p)) ;  
      colcurv = interp1([0; hicurve], u01, vmout4p) ;
    mout4p = reshape(colcurv,size(mout4p)) ;  
      colcurv = interp1([0; hicurve], u01, vmout1m) ;
    mout1m = reshape(colcurv,size(mout1m)) ;  
      colcurv = interp1([0; hicurve], u01, vmout2m) ;
    mout2m = reshape(colcurv,size(mout2m)) ;  
      colcurv = interp1([0; hicurve], u01, vmout3m) ;
    mout3m = reshape(colcurv,size(mout3m)) ;  
      colcurv = interp1([0; hicurve], u01, vmout4m) ;
    mout4m = reshape(colcurv,size(mout4m)) ;  
  end ;




  moutm = moutm .* mcrop ;
          %  Cut off part outside analysis circle
  mout1p = mout1p .* mcrop ;
  mout2p = mout2p .* mcrop ;
  mout3p = mout3p .* mcrop ;
  mout4p = mout4p .* mcrop ;
  mout1m = mout1m .* mcrop ;
  mout2m = mout2m .* mcrop ;
  mout3m = mout3m .* mcrop ;
  mout4m = mout4m .* mcrop ;


  figure(1) ;
  clf ;

  disp('Making Mean and Projections Plot') ;
  subplot(1,5,1) ;

  pcolor(mxgrid,mygrid,moutm) ;
    shading flat ;
    axis('square') ;
    title([imagstr surfstr ' mean']) ;
      ssstr = ['%SS,T=' num2str(100 * ssm / sst)] ;
    text(-2.4,3.6,ssstr) ;  
    if isurf == 2 | isurf == 3 ;   
          %  if heightmap - sphere or curvature
      caxis([0 1]) ;
    end ;
    colormap(comap) ;


  %  plot PC1 Components
  subplot(1,5,2) ;
  [kde, xgrid] = gpkde(pc1) ;
  plot(xgrid,kde,'w') ;
    title('Proj 1') ;
    vlim = axis ;
          %  Get axis values for drawing vertical lines
    left = min(xgrid) ;
    right = max(xgrid) ;
    range = right - left ;
    left = left - .05 * range ;
    right = right + .05 * range ;
    bottom = vlim(3) ;
    top = vlim(4) ;
    axis([left right bottom top]) ;
  hold on ;
    if igroup == 1 ;    %  Then color code bars for patients
      for ipat = 1:6 ;
          patcol = [(ipat > 3), (rem(ipat,4)> 1), rem(ipat,2)] ;
        set(gca,'ColorOrder',patcol) ;
              %  sets all colors to a default depending on ipat
              %  according to:  1 --->  0 0 1
              %                 2 --->  0 1 0
              %                 3 --->  0 1 1
              %                 4 --->  1 0 0
              %                 5 --->  1 0 1
              %                 6 --->  1 1 0
        ppc1 = pc1(vpatind == ipat) ;
              %  PC1's for current patient
        plot([ppc1 ppc1], [bottom top]) ;
      end ;
    else ;
      plot([pc1 pc1], [bottom top], 'c') ;
    end ;
    plot(xgrid,kde,'w') ;
          %  Do a replot to ensure seeing kde
    title('Proj, PC1') ;
  hold off ;


  %  plot PC2 Components
  subplot(1,5,3) ;
  [kde, xgrid] = gpkde(pc2) ;
  plot(xgrid,kde,'w') ;
    title('Proj 2') ;
    vlim = axis ;
          %  Get axis values for drawing vertical lines
    left = min(xgrid) ;
    right = max(xgrid) ;
    range = right - left ;
    left = left - .05 * range ;
    right = right + .05 * range ;
    bottom = vlim(3) ;
    top = vlim(4) ;
    axis([left right bottom top]) ;
  hold on ;
    if igroup == 1 ;    %  Then color code bars for patients
      for ipat = 1:6 ;
          patcol = [(ipat > 3), (rem(ipat,4)> 1), rem(ipat,2)] ;
        set(gca,'ColorOrder',patcol) ;
              %  sets all colors to a default depending on ipat
              %  according to:  1 --->  0 0 1
              %                 2 --->  0 1 0
              %                 3 --->  0 1 1
              %                 4 --->  1 0 0
              %                 5 --->  1 0 1
              %                 6 --->  1 1 0
        ppc2 = pc2(vpatind == ipat) ;
              %  PC2's for current patient
        plot([ppc2 ppc2], [bottom top]) ;
      end ;
    else ;
      plot([pc2 pc2], [bottom top], 'c') ;
    end ;
    plot(xgrid,kde,'w') ;
          %  Do a replot to ensure seeing kde
    title('Proj, PC2') ;
  hold off ;


  %  plot PC3 Components
  subplot(1,5,4) ;
  [kde, xgrid] = gpkde(pc3) ;
  plot(xgrid,kde,'w') ;
    title('Proj 3') ;
    vlim = axis ;
            %  Get axis values for drawing vertical lines
    left = min(xgrid) ;
    right = max(xgrid) ;
    range = right - left ;
    left = left - .05 * range ;
    right = right + .05 * range ;
    bottom = vlim(3) ;
    top = vlim(4) ;
    axis([left right bottom top]) ;
  hold on ;
    if igroup == 1 ;    %  Then color code bars for patients
      for ipat = 1:6 ;
          patcol = [(ipat > 3), (rem(ipat,4)> 1), rem(ipat,2)] ;
        set(gca,'ColorOrder',patcol) ;
              %  sets all colors to a default depending on ipat
              %  according to:  1 --->  0 0 1
              %                 2 --->  0 1 0
              %                 3 --->  0 1 1
              %                 4 --->  1 0 0
              %                 5 --->  1 0 1
              %                 6 --->  1 1 0
        ppc3 = pc3(vpatind == ipat) ;
              %  PC3's for current patient
        plot([ppc3 ppc3], [bottom top]) ;
      end ;
    else ;
      plot([pc3 pc3], [bottom top], 'c') ;
    end ;
    plot(xgrid,kde,'w') ;
            %  Do a replot to ensure seeing kde
    title('Proj, PC3') ;
  hold off ;
  

  %  plot PC4 Components
  subplot(1,5,5) ;
  [kde, xgrid] = gpkde(pc4) ;
  plot(xgrid,kde,'w') ;
    title('Proj 4') ;
    vlim = axis ;
          %  Get axis values for drawing vertical lines
    left = min(xgrid) ;
    right = max(xgrid) ;
    range = right - left ;
    left = left - .05 * range ;
    right = right + .05 * range ;
    bottom = vlim(3) ;
    top = vlim(4) ;
    axis([left right bottom top]) ;
  hold on ;
    if igroup == 1 ;    %  Then color code bars for patients
      for ipat = 1:6 ;
          patcol = [(ipat > 3), (rem(ipat,4)> 1), rem(ipat,2)] ;
        set(gca,'ColorOrder',patcol) ;
              %  sets all colors to a default depending on ipat
              %  according to:  1 --->  0 0 1
              %                 2 --->  0 1 0
              %                 3 --->  0 1 1
              %                 4 --->  1 0 0
              %                 5 --->  1 0 1
              %                 6 --->  1 1 0
        ppc4 = pc4(vpatind == ipat) ;
              %  PC4's for current patient
        plot([ppc4 ppc4], [bottom top]) ;
      end ;
    else ;
      plot([pc4 pc4], [bottom top], 'c') ;
    end ;
    plot(xgrid,kde,'w') ;
            %  Do a replot to ensure seeing kde
    title('Proj, PC4') ;
  hold off ;


    orient landscape ;
  eval(['print \matlab\steve\ps\' fsavstr '.ps -dpsc']) ;
          % make print of this first page


end ;





if imake == 2  | imake == 0 ;       %  Construct Movies 

  nframe = 21 ;
  nrpt = -5 ;

  disp('Making 1st eigenvector movie') ;
  figure(2) ;
  clf ;

  mmov1 = [] ;
  for iframe = 1:nframe ;
    w = (iframe - 1) / (nframe - 1) ;
    mout = (1 - w) * mout1p + w * mout1m ;

    colormap(comap) ;
    pcolor(mxgrid,mygrid,mout) ;
      shading flat ;
      axis('square') ;
          titstr = ['PC1'] ;
      title(titstr) ;
        if isurf == 2 | isurf == 3 ;   
            %  if heightmap - sphere or curvature
          caxis([0 1]) ;
        end ;

      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 ;

    mmov1(:,iframe) = getframe ;
  end ;

  movie(mmov1,nrpt) ;



  disp('Making 2nd eigenvector movie') ;
  figure(3) ;
  clf ;

  mmov2 = [] ;
  for iframe = 1:nframe ;
    w = (iframe - 1) / (nframe - 1) ;
    mout = (1 - w) * mout2p + w * mout2m ;

    colormap(comap) ;
    pcolor(mxgrid,mygrid,mout) ;
      shading flat ;
      axis('square') ;
          titstr = ['PC2'] ;
      title(titstr) ;
        if isurf == 2 | isurf == 3 ;   
            %  if heightmap - sphere or curvature
          caxis([0 1]) ;
        end ;

      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 ;

    mmov2(:,iframe) = getframe ;
  end ;

  movie(mmov2,nrpt) ;



  disp('Making 3rd eigenvector movie') ;
  figure(4) ;
  clf ;

  mmov3 = [] ;
  for iframe = 1:nframe ;
    w = (iframe - 1) / (nframe - 1) ;
    mout = (1 - w) * mout3p + w * mout3m ;

    colormap(comap) ;
    pcolor(mxgrid,mygrid,mout) ;
      shading flat ;
      axis('square') ;
          titstr = ['PC3'] ;
      title(titstr) ;
        if isurf == 2 | isurf == 3 ;   
            %  if heightmap - sphere or curvature
          caxis([0 1]) ;
        end ;

      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 ;

    mmov3(:,iframe) = getframe ;
  end ;

  movie(mmov3,nrpt) ;



  disp('Making 4th eigenvector movie') ;
  figure(5) ;
  clf ;

  mmov4 = [] ;
  for iframe = 1:nframe ;
    w = (iframe - 1) / (nframe - 1) ;
    mout = (1 - w) * mout4p + w * mout4m ;

    colormap(comap) ;
    pcolor(mxgrid,mygrid,mout) ;
      shading flat ;
      axis('square') ;
          titstr = ['PC4'] ;
      title(titstr) ;
        if isurf == 2 | isurf == 3 ;   
            %  if heightmap - sphere or curvature
          caxis([0 1]) ;
        end ;

      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 ;

    mmov4(:,iframe) = getframe ;
  end ;

  movie(mmov4,nrpt) ;


  eval(['save ' fsavstr ' mmov1 mmov2 mmov3 mmov4']) ;
          %  save in group indexed file

end ;




if imake == 3  | imake == 0 ;       %  Construct Control panel

  nrpt = -15 ;  
          %  number of repetitions (- gives back and forth)


  eval(['load ' fsavstr]) ;

  disp('Making Control Panel') ;
  figure(6) ;
  clf ;


  %  Create first push button: 1st eigenvector
  cbstr1 = '  figure(2) ;  colormap(comap) ; movie(mmov1,nrpt) ;' ;
  pushb1 = uicontrol(gcf,'Style','push','Position',[5 100 100 20],...
                'String','Princ. Comp. 1',...
                'CallBack',cbstr1) ;
                   %  "String" IS THE "LABEL" OF THE PUSH BUTTON

                   %  "CALLBACK" IS WHAT HAPPENS WHEN THE BUTTON IS 
                   %  PUSHED.  SINCE A LARGE NUMBER OF COMMANDS ARE
                   %  THEY ARE PUT INTO A SEPARATE SCRIPT FILE.

  %  Create second push button: 2nd eigenvector
  cbstr2 = '  figure(3) ;  colormap(comap) ; movie(mmov2,nrpt) ;' ;
  pushb2 = uicontrol(gcf,'Style','push','Position',[5 70 100 20],...
              'String','Princ. Comp. 2',...
              'CallBack',cbstr2) ;

  %  Create 3rd push button: 3rd eignevector
  cbstr3 = '  figure(4) ;  colormap(comap) ; movie(mmov3,nrpt) ;' ;
  pushb3 = uicontrol(gcf,'Style','push','Position',[5 40 100 20],...
              'String','Princ. Comp. 3',...
              'CallBack',cbstr3) ;

  %  Create 4th push button: 4th eigenvector
  cbstr4 = '  figure(5) ;  colormap(comap) ; movie(mmov4,nrpt) ;' ;
  pushb4 = uicontrol(gcf,'Style','push','Position',[5 10 100 20],...
              'String','Princ. Comp. 4',...
              'CallBack',cbstr4) ;


  disp('hit buttons for movies now') ;

end ;




if imake == 4  | imake == 0 ;       %  Convert files to mpeg

  eval(['load ' fsavstr]) ;

  disp('    turning first movie to mpeg') ;
    mpgwrite(mmov1,comap,['pc1' num2str(igroup) '.mpg']) ;
  disp('    turning second movie to mpeg') ;
    mpgwrite(mmov2,comap,['pc2' num2str(igroup) '.mpg']) ;
  disp('    turning third movie to mpeg') ;
    mpgwrite(mmov3,comap,['pc3' num2str(igroup) '.mpg']) ;
  disp('    turning fourth movie to mpeg') ;
    mpgwrite(mmov4,comap,['pc4' num2str(igroup) '.mpg']) ;
  disp('    done') ;


end ;
