disp('MATLAB Script File:  arvo1a.m') ;
%
%   For the Cornea Study, this makes color maps for ARVO poster
%         of April 1997
%   This is for means
%
format compact ;

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

ir0 = 2 ;     %  1 - 4.2
              %  2 - 1.7


ihires = 0 ;  %  Controls resolution of plot:
              %    0 - for lo res 64 x64  (fastest)
              %    1 - to interpolate to hi res 256 x 256  (slower)
              %    2 - for hi res 256 x 256  (REALLY slow)
              
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)


if ir0 == 1 ;
  r0 = 4.2 ;    %  Usual radius of analysis
elseif ir0 == 2 ;
  r0 = 1.7 ;    %  Usual radius of analysis
end ;


ngrid = 64 ; %  Construct surfaces as ngrid x ngrid array
if ihires == 2 ;
          %  then change ngrid to larger value
  ngrid = 256 ; %  Construct surfaces as ngrid x ngrid array
end ;
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 ;


for igp = 1:3 ;

  if igp == 1 ;
    igroup = 20 ;     %  20 - PRK Study (Pre Op)
  elseif igp == 2 ;
    igroup = 21 ;     %  21 - PRK Study (1 mo)
  elseif igp == 3 ;
    igroup = 212 ;     %  212 - PRK Study (12 mo)
  end ;



  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)
  %                                        (or "fm17")
  %    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 ;
    if ir0 == 1 ;
      eval(['load \matlab\steve\cornea\fm\fm42' fids '.mat']) ;
    elseif ir0 == 2 ;
      eval(['load \matlab\steve\cornea\fm\fm17' fids '.mat']) ;
    end ;
            %  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 ;





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



  d = size(mfv,1) ;
  n = size(mfv,2) ;
  %  get mean feature vector, and resids from that
  vmean = mean(mfv') ;


  fclose('all') ;


  %  Reshape representative guys to feature matrix form
  fmmean = reshape(vmean,19,7) ;


  surfstr = ' Curv''s' ;

  if ibas == 1 ;     %  then use the radial Fourier - Legendre basis
    [mhm,mhpm,mhppm] = corhpp(fmmean',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) ;
  moutm = moutm * 337.5 ;
          %  put on scale of diopters


  %  do explicit flip of x's and y's, since pcolor wants:
  %            rows as y's,   cols as x's
  moutm = moutm' ;


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

    vmoutm = moutm(:) ;
      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

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


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


  disp('Making Mean Plot') ;
  subplot(2,2,1) ;

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






  surfstr = ' Curv''s' ;

  fmmean = fmmean .* [ones(1,7);  zeros(18,7)] ;
          %  take only Fourier 0 terms

  if ibas == 1 ;     %  then use the radial Fourier - Legendre basis
    [mhm,mhpm,mhppm] = corhpp(fmmean',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) ;
  moutm = moutm * 337.5 ;
          %  put on scale of diopters

  %  do explicit flip of x's and y's, since pcolor wants:
  %            rows as y's,   cols as x's
  moutm = moutm' ;


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

    vmoutm = moutm(:) ;
      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

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

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


  disp('Making Mean Plot') ;
  subplot(2,2,3) ;

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



  if autop == 1 ;
    orient landscape ;
    eval(['print \matlab\steve\ps\arvo1a' num2str(ir0) ...
                                   num2str(igp) '.ps -dpsc']) ;
  end ;



end ;




