disp('MATLAB Script File:  corec12.m') ;
%
%   For the Cornea Study, this does simple reconstructions
%       from feature matrices, combining many on one plot
%       2nd Generation, using complete data sets
%   This modifies corec11.m, by using the feature matrices 
%       for the smaller diameter, r_0 = 1.7 (important for
%       the PRK data, instead of the usual r_0 = 4.2)
%
format compact ;

autop = 1 ;   %  1 to create a color postscript file

igroup = 212    %  1 - Repeatability Study (all)
                      %    Careful, this may not be done
              %  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)



r0 = 1.7 ;    %  Usual radius of analysis

ngrid = 64 ; %  Construct surfaces as ngrid x ngrid array
hrngrid = 256 ; %  Interpolate to hrngrid x hrngrid for 
                %  higher resolution part

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 ' ;
  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 "fm17" 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\fm17' 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 "fm17" 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\fm17' 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']) ;




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


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

npage = ceil(n / 12) ;
          %  number of pages needed to print all of these


for ipage = 1:npage ;

  if ipage < npage ;
    npic = 12 ;
  else ;
    npic = n - 12 * (ipage - 1) ;
  end ;

  for ipic = 1:npic ;

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

    iind = ipic + 12 * (ipage - 1) ;
    fv = mfv(:,iind) ;


    %  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(3,4,ipic) ;

    pcolor(mxgrid,mygrid,mout) ;
      shading flat ;
      axis('square') ;
        if igroup == 1 ;
          titstr = [imagstr ' pat' num2str(vpatind(iind)) ...
                            '  ' num2str(veyeind(iind)) ] ;
        elseif igroup < 300 ;
          titstr = [imagstr ' pat' num2str(vpatind(iind)) ...
                            '  ' num2str(veyeind(iind)) ] ;
        else ;
          titstr = [imagstr ' pat' num2str(vpatind(iind)) ...
                            '  ' num2str(veyeind(iind)) ] ;
       end ;
      title(titstr) ;
        if isurf == 2 | isurf == 3 ;   
          %  if heightmap - sphere or curvature
          caxis([0 1]) ;
        end ;

  end ;


  if autop == 1 ;
    orient landscape ;
    eval(['print \matlab\steve\ps\co' num2str(igroup) ...
                    num2str(ipage) '.ps -dpsc']) ;
  end ;

end ;



