disp('MATLAB Script File:  cofmap1.m') ;
%
%   For the Cornea Study, this makes gray scale maps of
%                     feature matrices
%   Copied from corecon6.m
%
format compact ;

autop = 1 ;   %  1 to create a postscript file


iimag = 125 ;  
              %  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


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 ;
  %  Convert c-s to 10 x 7 amplitude - phase form, 
  %        using block from cofmtst1.m
  mamp = abs(fm(1,:)) ;
  mpha = atan2(0,fm(1,:)) * 180 / pi ;
  for i = 1:9 ;
      vamp = sqrt(fm(2*i,:).^2 + fm(2*i+1,:).^2) ;
    mamp = [mamp; vamp] ;
      vpha = atan2(fm(2*i+1,:),fm(2*i,:)) * 180 / pi ;
        vind = find(vpha < 0) ;
      vpha(vind) = vpha(vind) + 360 ;
    mpha = [mpha; vpha] ;
  end ;

  %  Feature matrix form for output
  fmamp = fm(1,:) ;   %  To get signs right
  fmamp = [fmamp; mamp(2:10,:)] ;  %  Use amplitude for rest
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 ;

  %  Feature matrix form for output
  fmamp = mamp(1,:) .* cos(mpha(1,:) * pi / 180) ;
          %  Here Jinting used names mnamp & mnpha
  fmamp = [fmamp; mamp(2:10,:)] ;  %  Use amplitude for rest
end ;

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


%  Look only at magnitudes
fmamp = abs(fmamp) ;

%  Put on scale 0-64 for image construction
mout = 64 * fmamp / max(max(fmamp)) ;

subplot(2,1,1) ;
  image(0:6,0:9,mout) ;
    colormap('gray') ;
    colorbar ;
      titstr = ['Abs(FM), ' imagstr] ;
    title(titstr) ;

%  Look at logs
mout = log10(fmamp + 10^(-20)) ;
          %  Avoids problems with 0's in log
mout = mout + 4 ;
          %  Puts 0 at -4 (thus truncates stuff below 10^(-4) in image)
mout = 64 * mout / max(max(mout)) ;

subplot(2,1,2) ;
  image(0:6,0:9,mout) ;
    colormap('gray') ;
    colorbar ;
      titstr = ['log(Abs(FM)>.0001), ' imagstr] ;
    title(titstr) ;


  orient tall ;
if autop == 1 ;
  eval(['print \matlab\steve\ps\fmap' num2str(iimag) '.ps -dps']) ;
end ;

