disp('MATLAB Script File:  corecon6.m') ;
%
%   For the Cornea Study, this reconstructs surfaces from
%                     feature matrices
%       Modification of corecon3.m, and corecon5.m
%       which it improves by:
%            (i) reading from .mat files,
%            (ii) solving orientation problems
%            (iii) combining individs, and means into one program
%
format compact ;

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

ihires = 2 ;  %  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)

iimag = 107 ;  
              %  1 - Individual 84_1
              %  2 - Individual 85_1
              %  3 - Individual 122_2
              %  4 - Individual 121_2
              %  5 - Individual 55_1
              %  6 - Individual 28_2
              %  7 - Individual 130_1
              %  8 - Individual 139_2
              %  9 - Individual 6_1
              %  10 - Individual 53_1
              %  11 - Individual 43_1

              %  101 - Mean Pre
              %  102 - Mean Pre Right
              %  103 - Mean Pre Left
              %  104 - Mean Diag Ker
              %  105 - Mean Diag Ker Right
              %  106 - Mean Diag Ker Left
              %  107 - Mean Ksuspect
              %  108 - Mean Ksuspect Right
              %  109 - Mean Ksuspect Left
              %  110 - Mean 1 Month
              %  111 - Mean 1 Month Right
              %  112 - Mean 1 Month Left
              %  113 - Mean 6 Month
              %  114 - Mean 6 Month Right
              %  115 - Mean 6 Month Left
              %  116 - Mean 1 Year
              %  117 - Mean 1 Year Right
              %  118 - Mean 1 Year Left
              %  119 - Mean All K
              %  120 - Mean All K Right
              %  121 - Mean All K Left
              %  122 - Mean all P
              %  123 - Mean all P Right
              %  124 - Mean all P Left
              %  125 - Mean Overall
              %  126 - Mean Overall Right
              %  127 - Mean Overall Left

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)

%figure(1) ;
%clg reset ;
%disp('Hit any key to continue after clearing graphics') ;
%pause ;

r0 = 4.2 ;    %  Usual radius of analysis
%spherad = 9.3 ;%  Radius of sphere to subtract
%  now set individually for each image

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 ;

if iimag == 1 ;
  imagstr = 'Image 84_1' ;
  fmistr = 'fm420841.mat' ;
  pdist = .123 ;
  pang = 123 ;
  spherad = 7.941 ;  %  Radius of sphere to subtract
  cstep = 6 ;  %  Step size for color scale
elseif iimag == 2 ;
  imagstr = 'Image 85_1' ;
  fmistr = 'fm420851.mat' ;
  pdist = .254 ;
  pang = 68 ;
  spherad = 7.258 ; %  Radius of sphere to subtract
  cstep = 25 ;  %  Step size for color scale
elseif iimag == 3 ;
  imagstr = 'Image 122_2' ;
  fmistr = 'fm421222.mat' ;
  pdist = .686 ;
  pang = 89 ;
  spherad = 7.8488 ; %  Radius of sphere to subtract
  cstep = 35 ;  %  Step size for color scale
elseif iimag == 4 ;
  imagstr = 'Image 121_2' ;
  fmistr = 'fm421212.mat' ;
  pdist = .525 ;
  pang = 108 ;
  spherad = 337.5 / 41 ;  %  Radius of sphere to subtract (Diotpers)
  cstep = 40 ;  %  Step size for color scale
elseif iimag == 5 ;
  imagstr = 'Image 55_1' ;
  fmistr = 'fm420551.mat' ;
  pdist = .439 ;
  pang = 121 ;
  spherad = 337.5 / 43.5 ;  %  Radius of sphere to subtract (Diotpers)
  cstep = 12 ;  %  Step size for color scale
elseif iimag == 6 ;
  imagstr = 'Image 28_2' ;
  fmistr = 'fm420282.mat' ;
  pdist = .426 ;
  pang = 115 ;
  spherad = 8.0357 ; %  Radius of sphere to subtract
  cstep = 5 ;  %  Step size for color scale
elseif iimag == 7 ;
  imagstr = 'Image 130_1' ;
  fmistr = 'fm421301.mat' ;
  pdist = .458 ;
  pang = 112 ;
  spherad = 337.5 / 39.7 ;  %  Radius of sphere to subtract (Diotpers)
  cstep = 10 ;  %  Step size for color scale
%  spherad = 337.5 / 37.9 ;  %  Radius of sphere to subtract (Diotpers)
elseif iimag == 8 ;
  imagstr = 'Image 139_2' ;
  fmistr = 'fm421392.mat' ;
  pdist = .417 ;
  pang = 114 ;
  spherad = 8.4375 ; %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 9 ;
  imagstr = 'Image 6_1' ;
  fmistr = 'fm420061.mat' ;
  pdist = .454 ;
  pang = 116 ;
  spherad = 7.7586 ; %  Radius of sphere to subtract
  cstep = 6 ;  %  Step size for color scale
elseif iimag == 10 ;
  imagstr = 'Image 53_1' ;
  fmistr = 'fm420531.mat' ;
  pdist = .375 ;
  pang = 113 ;
  spherad = 8.766 ; %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 11 ;
  imagstr = 'Image 43_1' ;
  fmistr = 'fm420431.mat' ;
  pdist = .557 ;
  pang = 91 ;
  spherad = 8.5443 ; %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale


elseif iimag == 101 ;
  imagstr = 'Mean Pre' ;
  fmistr = 'mn42pre.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 102 ;
  imagstr = 'Mean Pre Right' ;
  fmistr = 'mn42pre.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 103 ;
  imagstr = 'Mean Pre Left' ;
  fmistr = 'mn42pre.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 104 ;
  imagstr = 'Mean Diag Ker' ;
          %  Diagnosed Keratoconus 
          %          - called ker by NKT
          %          - called kdi by JTZ
  fmistr = 'mn42kdi.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 105 ;
  imagstr = 'Mean Diag Ker Right' ;
  fmistr = 'mn42kdi.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 106 ;
  imagstr = 'Mean Diag Ker Left' ;
  fmistr = 'mn42kdi.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 107 ;
  imagstr = 'Mean Ksuspect' ;
  fmistr = 'mn42ksu.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 108 ;
  imagstr = 'Mean Ksuspect Right' ;
  fmistr = 'mn42ksu.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 109 ;
  imagstr = 'Mean Ksuspect Left' ;
  fmistr = 'mn42ksu.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 110 ;
  imagstr = 'Mean 1 Month' ;
  fmistr = 'mn421mo.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 111 ;
  imagstr = 'Mean 1 Month Right' ;
  fmistr = 'mn421mo.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 112 ;
  imagstr = 'Mean 1 Month Left' ;
  fmistr = 'mn421mo.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 113 ;
  imagstr = 'Mean 6 Month' ;
  fmistr = 'mn426mo.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 114 ;
  imagstr = 'Mean 6 Month Right' ;
  fmistr = 'mn426mo.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 115 ;
  imagstr = 'Mean 6 Month Left' ;
  fmistr = 'mn426mo.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 116 ;
  imagstr = 'Mean 1 Year' ;
  fmistr = 'mn421yr.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 117 ;
  imagstr = 'Mean 1 Year Right' ;
  fmistr = 'mn421yr.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 118 ;
  imagstr = 'Mean 1 Year Left' ;
  fmistr = 'mn421yr.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 119 ;
  imagstr = 'Mean All K Study' ;
          %  All 3 groups in Keratoconus study
          %          - called kstud by NKT
          %          - called ker by JTZ
  fmistr = 'mn42ker.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 120 ;
  imagstr = 'Mean All K Study Right' ;
  fmistr = 'mn42ker.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 121 ;
  imagstr = 'Mean All K Study Left' ;
  fmistr = 'mn42ker.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 122 ;
  imagstr = 'Mean All P Study' ;
          %  All 4 groups in PRK study
          %          - called pstud by NKT
          %          - called pos by JTZ
          %                   nor by JTZ are not used (same as pre).
  fmistr = 'mn42pos.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 123 ;
  imagstr = 'Mean All P Study Right' ;
  fmistr = 'mn42pos.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 124 ;
  imagstr = 'Mean All P Study Left' ;
  fmistr = 'mn42pos.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 125 ;
  imagstr = 'Overall Mean' ;
          %  Note:  pre group is included twice as "pre" and as "nor"
  fmistr = 'mn42all.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 126 ;
  imagstr = 'Overall Mean Right' ;
  fmistr = 'mn42all.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 127 ;
  imagstr = 'Overall Mean Left' ;
  fmistr = 'mn42all.mat' ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
end ;

  
['load \matlab\steve\cornea\fm\' fmistr]
eval(['load \matlab\steve\cornea\fm\' fmistr]) ;



if iimag < 100 ;
  mfeat = fm' ;    %  Jinting used name fm
else ;
  if rem(iimag,3) == 2 ;   %  Then this is an "all" image
    mamp = mnamp ;
    mpha = mnpha ;
  elseif rem(iimag,3) == 0 ;   %  Then this is a "right" image
    mamp = mnampr ;
    mpha = mnphar ;
  elseif rem(iimag,3) == 1 ;   %  Then this is a "left" image
    mamp = mnampl ;
    mpha = mnphal ;
  end ;

  fm = mamp(1,:) .* cos(mpha(1,:) * pi / 180) ;
          %  Here Jinting used names mnamp & mnpha
  for i = 2:10 ;
    fm = [fm; mamp(i,:) .* cos(mpha(i,:) * pi / 180)] ;
    fm = [fm; mamp(i,:) .* sin(mpha(i,:) * pi / 180)] ;
  end ;
  mfeat = fm' ;
end ;


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




disp('Generating surface') ; 
if isurf == 1 | isurf == 2 ;   %  then work with height map
  surfstr = ' Heights' ;

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

    mout = corh(mfeat,mrgrid,mthetag,r0) ;
          %  Radial Fourier - Legendre reconstruction 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 = ' Curvatures' ;

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

    [mh,mhp,mhpp] = corhpp(mfeat,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 ;




%  mout = hh + sphht ;
           %  REPLACES WITH JINTING'S HEIGHTS, FOR TESTING, HIS
           %  TESTHH STUFF IS IMAGE 28-2.




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



if ihires == 1 ;
          %  Then compute higher resolution stuff
  hrxg = linspace(-r0,r0,hrngrid) ;
  hryg = hrxg ;
  [mhrxg,mhryg] = meshgrid(hrxg,hryg) ;
  mhrout = interp2(mxgrid,mygrid,mout,mhrxg,mhryg) ;
  mcrop = ones(hrngrid) ;            
  vind = find(mhrxg.^2 + mhryg.^2 > r0^2) ;
          %  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.       
  mhrout = mhrout .* mcrop ;
end ;






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

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

    vmout = mout(:) ;
          %  vector version for easier manipulation
          %  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 ;




if ihires == 1 ;
          %  Then compute higher resolution stuff
  hrxg = linspace(-r0,r0,hrngrid) ;
  hryg = hrxg ;
  [mhrxg,mhryg] = meshgrid(hrxg,hryg) ;
  mhrout = interp2(mxgrid,mygrid,mout,mhrxg,mhryg) ;
  mcrop = ones(hrngrid) ;            
  vind = find(mhrxg.^2 + mhryg.^2 > r0^2) ;
          %  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.       
  mhrout = mhrout .* mcrop ;
else ;
  mout = mout .* mcrop ;
          %  Cut off part outside analysis circle
end ;



disp('Making Plot') ;

if ihires == 1 ;
  pcolor(mhrxg,mhryg,mhrout) ;
else ;
  pcolor(mxgrid,mygrid,mout) ;
end ;
  shading flat ;
  axis('square') ;
  title([imagstr ' Recon' surfstr]) ;
  if isurf == 2 | isurf == 3 ;   
          %  if heightmap - sphere or curvature
    caxis([0 1]) ;
    if iimag < 100 ;  %  Then have individual, so add Keratron Center
      hold on ;
        pangr = (2 * pi * pang) / 360 ;
          %  Convert to radians
        polar((pi + pangr),pdist,'*w') ;
            %  Show Keratron center
        axis('off') ;
      hold off ;
%  careful about these for left eye (because feature matrices 
%      have been "flipped").
    end ;
  end ;


if autop == 1 ;
  if isurf == 2 ;
    eval(['print \matlab\steve\ps\re6' num2str(iimag) 'h.ps -dpsc']) ;
  elseif isurf == 3 ;
    eval(['print \matlab\steve\ps\re6' num2str(iimag) 'c.ps -dpsc']) ;
  end ;
end ;

