disp('MATLAB Script File:  comov1.m') ;
%
%   For the Cornea Study, this does a movie 
%       morph'ing through time for the PRK data
%       2nd Generation, using complete data sets

format compact ;

imake = 0 ;      %  0 - do everything
                 %  1 - make figure 1, static plots
                 %  2 - make figure 2, with movies
                 %  3 - make figure 3 a GUI menu, and run movies in fig. 2
                 %  4 - make two mpeg files

ipat = 109 ;      %  can be 100,...,109, patients in PRK data
                  %  with all 4 images present
                  %    most interesting to me is: 109
                  %  Nancy wants: 100 (not so interesting),
                  %               107 (very good),
                  %               108 (very good, esp. when comp. to 107)

%igroup = 212    %  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)

              
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)




%  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) ;



%  Set up Figure windows
%
fhand1 = figure(1) ;
  set(fhand1,'Position',[330,200,300,300]) ;

fhand2 = figure(2) ;
  set(fhand2,'Position',[175,350,150,150]) ;

fhand3 = figure(3) ;
  set(fhand3,'Position',[400,20,220,90]) ;



%  Set file save string
%
fsavstr = ['cm1p' num2str(ipat)] ;


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 ;



            %  following block of code from corec11.m
  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



  ipatstr = num2str(ipat) ;
  patstr = ['P ' ipatstr] ;
    disp(['Loading feature matrices for ' patstr]) ;
    vflag = (vpatind == ipat) ;

  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 ;




  fclose('all') ;


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





  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 \matlab\steve\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


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

  end ;




  d = size(mfv,1) ;
  n = size(mfv,2) ;



  %  make figure of all 4 next to each other
  figure(1) ;
  for ipic = 1:4 ;

    disp(['  working on pic ' num2str(ipic)]) ;

    fv = mfv(:,ipic) ;


    %  Reshape representative guys to feature matrix form
    fm = reshape(fv,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

        mout = corh(fm,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


        mout = mout + 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

        [mh,mhp,mhpp] = corhpp(fm',mrgrid,mthetag,r0) ;

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

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

      mout = -mhpp ./ (1 + mhp.^2).^(3/2) ;
      mout = mout * 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
    mout = mout' ;





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

      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.
      mout = reshape(colcurv,size(mout)) ;  


    elseif isurf == 3 ;

      vmout = mout(:) ;
            %  NOTE: interp1 seems to need vector argument
        indlo = find(vmout < 0) ;
      vmout(indlo) = 0 * ones(length(indlo),1) ;
        indup = find(vmout > max(hicurve)) ;
      vmout(indup) = max(hicurve) * ones(length(indup),1) ;
            %  NOTE: interp1 seems to need values within grid

      colormap(comap) ;

        u01 = linspace(0,1,ncomap+1) ;
        colcurv = interp1([0; hicurve], 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.
          %  NOTE: interp1 seems to need vector argument
      mout = reshape(colcurv,size(mout)) ;  
    end ;




    mout = mout .* mcrop ;
          %  Cut off part outside analysis circle



    subplot(2,2,ipic) ;

    pcolor(mxgrid,mygrid,mout) ;
      shading flat ;
      axis('square') ;
          titstr = ['PRK ' patstr ',  ' num2str(vtime(ipic)) ' mo.'] ;
      title(titstr) ;
        if isurf == 2 | isurf == 3 ;   
            %  if heightmap - sphere or curvature
          caxis([0 1]) ;
        end ;

    if ipic == 2 ;
      mout2 = mout ;
    elseif ipic == 3 ;
      mout3 = mout ;
    elseif ipic == 4 ;
      mout4 = mout ;
    end ;

  end ;


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

end ;




%  Now work on movie construction
%
if imake == 2  | imake == 0 ;
  disp('making first movie') ;

  figure(2) ;
    clf ;

%  if imake == 0 ;
%    disp('adjust these window sizes, and hit any key') ;
%    pause ;
%  end ;
  nframe = 21 ;
  nrpt = -5 ;

  mmov1 = [] ;
%  mmov1 = moviein(nframe) ;
  mmov3 = [] ;
  for iframe = 1:nframe ;
    w = (iframe - 1) / (nframe - 1) ;
    mout = (1 - w) * mout2 + w * mout3 ;

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

      hold on ;   % add bar showing time gone
        bars = 2.5 ;
        barl = 1.5 ;
        bary = 4 ;
        plot([bars,bars+w*barl],[bary,bary],'-y') ;
          vhand = get(gca,'Children') ;
          set(vhand(1),'LineWidth',5) ;
      hold off ;

      if iframe == 1 ;
        text(-3.95,3.7,'1 Mo.') ;
      elseif  iframe == nframe ;
        text(-3.95,3.7,'6 Mo.') ;
      end ;    

    mmov1(:,iframe) = getframe ;
    mmov3(:,iframe) = getframe ;
  end ;

  movie(mmov1,nrpt) ;




  disp('making 2nd movie') ;

  figure(2) ;

  mmov2 = [] ;
%  mmov2 = moviein(nframe) ;
  for iframe = 1:nframe ;
    w = (iframe - 1) / (nframe - 1) ;
    mout = (1 - w) * mout3 + w * mout4 ;

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

      hold on ;   % add bar showing time gone
        bars = 2.5 ;
        barl = 1.5 ;
        bary = 4 ;
        plot([bars,bars+w*barl],[bary,bary],'-c') ;
          vhand = get(gca,'Children') ;
          set(vhand(1),'LineWidth',5) ;
      hold off ;

      if iframe == 1 ;
        text(-3.95,3.7,'6 Mo.') ;
      elseif  iframe == nframe ;
        text(-3.95,3.7,'12 Mo.') ;
      end ;    

    mmov2(:,iframe) = getframe ;
    mmov3(:,nframe+iframe) = getframe ;
  end ;

  movie(mmov2,nrpt) ;


  disp('  hit any key to see the joint movie') ;
  pause ;

  movie(mmov3,nrpt) ;


  eval(['save ' fsavstr ' mmov1 mmov2 mmov3']) ;


end ;




%  Now work on movie showing
%
if imake == 3  | imake == 0 ;

  nrpt = -15 ;  
          %  number of repetitions (- gives back and forth)
  fps = 4 ;
          %  frames per second (default = 12)

  eval(['load ' fsavstr]) ;

  disp('building user interface') ;

  figure(3) ;
    clf ;

  %  Create frame
  %
  mainh = uicontrol(gcf,'Style','frame',...
               'Position',[5 2 200 80]) ; 
                     %  NUMBER OF PIXELS FROM LEFT EDGE
                     %  NUMBER OF PIXELS FROM BOTTOM
                     %  NUMBER OF PIXELS IN WIDTH
                     %  NUMBER OF PIXELS IN HEIGHT

%  %  Create first push button: 1 morph of 1-6 mo.
%  cbstr1 = '  figure(2) ;  colormap(comap) ; movie(mmov1,1,fps) ;' ;
%  pushb1 = uicontrol(gcf,'Style','push','Position',[10 80 190 20],...
%              'String','Slow Single Morph, 1 - 6 mo.',...
%              '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: 1 morph of 6-12 mo.
%  cbstr2 = '  figure(2) ;  colormap(comap) ; movie(mmov2,1,fps) ;' ;
%  pushb2 = uicontrol(gcf,'Style','push','Position',[10 55 190 20],...
%              'String','Slow Single Morph, 6 - 12 mo.',...
%              'CallBack',cbstr2) ;


  %  Create 1st push button: repeated morphs of 1-6 mo.
  cbstr1 = '  figure(2) ;  colormap(comap) ; movie(mmov1,nrpt) ;' ;
  pushb1 = uicontrol(gcf,'Style','push','Position',[10 55 190 20],...
              'String','Repeated Morphs, 1 - 6 mo.',...
              'CallBack',cbstr1) ;

  %  Create 2nd push button: repeated morphs of 6-12 mo.
  cbstr2 = '  figure(2) ;  colormap(comap) ; movie(mmov2,nrpt) ;' ;
  pushb2 = uicontrol(gcf,'Style','push','Position',[10 30 190 20],...
              'String','Repeated Morphs, 6 - 12 mo.',...
              'CallBack',cbstr2) ;

  %  Create 3rd push button: repeated morphs of 1 - 12 m0.
  cbstr3 = '  figure(2) ;  colormap(comap) ; movie(mmov3,nrpt) ;' ;
  pushb3 = uicontrol(gcf,'Style','push','Position',[10 5 190 20],...
              'String','Repeated Morphs, 1 - 12 mo.',...
              'CallBack',cbstr3) ;



  disp('hit buttons for movies now') ;


end ;



%  Now work on saving movie as MPEG
%
if imake == 4  | imake == 0 ;

  eval(['load ' fsavstr]) ;

  disp('    turning first movie to mpeg') ;
    mpgwrite(mmov1,comap,['cm1' num2str(ipat) '.mpg']) ;
  disp('    turning second movie to mpeg') ;
    mpgwrite(mmov2,comap,['cm2' num2str(ipat) '.mpg']) ;
  disp('    turning third movie to mpeg') ;
    mpgwrite(mmov3,comap,['cm3' num2str(ipat) '.mpg']) ;
  disp('    done') ;


end ;


