disp('MATLAB Script File:  cocd3.m') ;
%
%   For the Cornea Study, this does a "curves as data" approach
%       to understanding the variability in the corneal groups.
%   Copied from cocd3.m, which it modifies by looking at variability
%       about a "central image"
%
format compact ;

autop = 1 ;   %  1 to create a postscript file

igroup = 1 ;  %  1 - 1 Month      (Central Image = "All P")
              %  2 - Kerataconic  (Central Image = "All K")

idisp = 1 ;   %  1 - Fourier 0 and 1
              %  2 - Fourier 0 only
              %  3 - Fourier 1 only
              %  4 - Fourier 2 only
              %  5 - Fourier 3 only
              %  6 - Fourier 4 only
              %  7 - 3 x 3 block for Ker
              %  8 - 5 x 2 block for Ker

isurf = 3 ;   %  1 - for height map
              %  2 - for height map - sphere  (not working)
              %  3 - for curvature map

ibas = 1 ;    %  1 - Radial Fourier - Legendre
              %  2 - Diameter Fourier - Legendre
              %            (not working yet)
              %  3 - Zernike
              %            (not working yet)


if idisp == 1 ;   %  Fourier 0 & 1
  mind = [     0,0; 0,1; 0,2;
          0,3; 0,4; 0,5; 0,6;
          1,1; 1,2; 1,3      ] ;
          %  This needs to be a 10 x 2 matrix, with entries
          %  indicating features to show (Fourier/Lengedre)
          %  Note upper left is "everything" and 
          %      lower right is "final residual".
          %  Entries are placed to show where plot will appear.
elseif idisp == 2 ;   %  Fourier 0 only
  mind = [     0,0; 0,1; 0,2;
          0,3; 0,4; 0,5; 0,6
                             ] ;
elseif idisp == 3 ;   %  Fourier 1 only
  mind = [     1,0; 1,1; 1,2;
          1,3; 1,4; 1,5; 1,6
                             ] ;
elseif idisp == 4 ;   %  Fourier 2 only
  mind = [     2,0; 2,1; 2,2;
          2,3; 2,4; 2,5; 2,6
                             ] ;
elseif idisp == 5 ;   %  Fourier 3 only
  mind = [     3,0; 3,1; 3,2;
          3,3; 3,4; 3,5; 3,6
                             ] ;
elseif idisp == 6 ;   %  Fourier 4 only
  mind = [     4,0; 4,1; 4,2;
          4,3; 4,4; 4,5; 4,6
                             ] ;
elseif idisp == 7 ;   %  3 x 3 block for Ker
  mind = [     0,0; 0,1; 0,2;
          1,0; 1,1; 1,2; 2,0;
          2,1; 2,2;          ] ;
elseif idisp == 8 ;   %  5 x 2 block for Ker
  mind = [     0,0; 0,1; 1,0;
          1,1; 2,0; 2,1; 3,1;
          3,2; 4,1; 4,2      ] ;
end ;
vindoa = 10 * mind(:,2) + mind(:,1) + 1 ;

r0 = 4.2 ;    %  Usual radius of analysis

ngrid = 100 ; 

nr = 10 ;
nc = 7 ;

if igroup == 1 ;
  namestr = '\matlab\steve\cornea\fm\fname1mo.asc' ;
  imagstr = '1 Mo.' ;
  icent = 122 ;    %  cofmap1.m index for Centering Feature Matrix
              %  122 - Mean all P
  centstr = 'all P' ;
elseif igroup == 2 ;
  namestr = '\matlab\steve\cornea\fm\fnameker.asc' ;
  imagstr = 'Ker.' ;
  icent = 119 ;    %  cofmap1.m index for Centering Feature Matrix
              %  119 - Mean All K
  centstr = 'all K' ;
end ;

disp('Loading feature matrices for this group') ;
fid = fopen(namestr,'r') ;
vfindx = fscanf(fid,'%4d') ;

prestr = '\matlab\steve\cornea\fm\fm42' ;
mfa = [] ;
for i = 1:size(vfindx,1) ;

  if vfindx(i) < 10 ;
    eval(['load ' prestr '000' num2str(vfindx(i))]) ;
  elseif vfindx(i) < 100 ;
    eval(['load ' prestr '00' num2str(vfindx(i))]) ;
  elseif vfindx(i) < 1000 ;
    eval(['load ' prestr '0' num2str(vfindx(i))]) ;
  else ;
    eval(['load ' prestr num2str(vfindx(i))]) ;
  end ;
          %  This creates a number of variables, especially fm,
          %  a 19 x 7 "feature matrix"

  %  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
  mfa = [mfa fmamp(:)] ;
          %  put matrix into a column vector and add to right
end ;


disp('Load centering feature matrix') ;
if icent == 119 ;
%  imagstr = 'Mean All K Study' ;
          %  All 3 groups in Keratoconus study
          %          - called kstud by NKT
          %          - called ker by JTZ
  fmistr = 'mn42ker.mat' ;
elseif icent == 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' ;
end ;

['load \matlab\steve\cornea\fm\' fmistr]
eval(['load \matlab\steve\cornea\fm\' fmistr]) ;
mamp = mnamp ;
mpha = mnpha ;
%  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

mfa = mfa - vec2mat(fmamp(:),size(vfindx,1)) ;


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

disp(['feature matrix has size ' num2str(d) ' rows ']) ;
disp(['and ' num2str(n) ' cols']) ;


fclose('all') ;


rgrid = linspace(0,r0,ngrid)' ;
          %  Column vector of radii

%  Generate matrix of Legendre's
ell0 = 6 ;
mpolyarg = (2 * rgrid ./ r0) - 1 ;

%  Do ell = 0 Legendre
llleg = ones(size(mpolyarg)) ;
          %  Last, last Legendre for future updating
lllegp = zeros(size(mpolyarg)) ;
          %  Last, last Legendre prime for future updating
lllegpp = zeros(size(mpolyarg)) ;
          %  Last, last Legendre prime, prime for future updating
mleg = llleg ;
mlegp = lllegp ;
mlegpp = lllegpp ;

%  Do ell = 1 Legendre
lleg = mpolyarg ;
          %  Last Legendre for future updating
llegp = ones(size(mpolyarg)) ;
          %  Last Legendre prime for future updating
llegpp = zeros(size(mpolyarg)) ;
          %  Last Legendre prime, prime for future updating
mleg = [mleg lleg] ;
mlegp = [mlegp llegp] ;
mlegpp = [mlegpp llegpp] ;

%  Loop through rest of Legendre's
for ell = 2:ell0 ;
  leg = (1/ell) .* ((2*ell-1) .* mpolyarg .* lleg - (ell-1) .* llleg) ;
              %  Legendre update formula
  legp = (1/ell) .* ((2*ell-1) .* (lleg + mpolyarg .* llegp) ...
                            - (ell-1) .* lllegp) ;
          %  Legendre prime update formula
  legpp = (1/ell) .* ((2*ell-1) .* (2 * llegp + mpolyarg .* llegpp) ...
                            - (ell-1) .* lllegpp) ;
          %  Legendre prime, prime update formula
  mleg = [mleg leg] ;
  mlegp = [mlegp legp] ;
  mlegpp = [mlegpp legpp] ;

  llleg = lleg ;
  lleg = leg ;
  lllegp = llegp ;
  llegp = legp ;
  lllegpp = llegpp ;
  llegpp = legpp ;
          %  Update for next run through loop
end ;

mlegp = mlegp * (2 ./ r0) ;
mlegpp = mlegpp * (2 ./ r0).^2 ;



disp('Reconstructing Individual Amplitudes') ;
rdat = [] ;
for i = 1:n ;    %  Loop through individuals
  vlegcoef = mfa(:,i) ;

  vdat = zeros(ngrid,1) ;
  for ifour = 1:10 ;    %  Loop through Fourier levels
    legcoef = vlegcoef((0:10:60)+ifour) ;

    %  reconstruct Legendre's 
    mrecon = mleg' .* vec2mat(legcoef,ngrid) ;
    mrecon = sum(mrecon)' ;
    mreconp = mlegp' .* vec2mat(legcoef,ngrid) ;
    mreconp = sum(mreconp)' ;
    mreconpp = mlegpp' .* vec2mat(legcoef,ngrid) ;
    mreconpp = sum(mreconpp)' ;
  
    mcurv = -mreconpp ./ (1 + mreconp.^2).^(3/2) ;
    mcurv = mcurv * 337.5 ;
          %  put on scale of diopters

    if isurf == 1 ;
      vdat = vdat + mrecon ;
    elseif isurf == 3 ;
      vdat = vdat + mcurv ;
    end ;
  end ;

  rdat = [rdat vdat] ;
end ;


sst = sum(sum(mfa.^2)) ;

%  Reset Graphics to get proper blanks
clf ;

disp('   Making 1st plot') ;
textx = 0.2 ;
if isurf == 1 ;
  bottom = -1.2 ;
  top = .3 ;
elseif isurf == 3 ;
  bottom = -40 ;
  top = 70 ;
end ;
subplot(3,4,1) ;
  plot(rgrid,rdat) ;
    title([imagstr ', Raw - ' centstr]) ;
    axis([0 r0 bottom top]) ;
      ssstr = ['SS,T = ' num2str(sst)] ;
      texty = bottom + .9 * (top - bottom) ;
    text(textx,texty,ssstr) ;  

  
%  Loop through single direction plots
ssr = sst ;
rmresid = rdat ;
for ip = 2:(length(vindoa)+1) ;
  disp(['    Making plot' num2str(ip)]) ;
  indoa = vindoa(ip - 1) ;
  legcoef = mfa(indoa,:) ;
  spread = std(legcoef) ;
  ssoa = sum((legcoef.^2)') ;
  ssr = ssr - ssoa ;

  %  reconstruct Legendre's 
  mrecon = vec2mat(mleg(:,mind(ip-1,2)+1),n) .* ...
                                     vec2mat(legcoef,ngrid) ;
  mreconp = vec2mat(mlegp(:,mind(ip-1,2)+1),n) .* ...
                                     vec2mat(legcoef,ngrid) ;
  mreconpp = vec2mat(mlegpp(:,mind(ip-1,2)+1),n) .* ...
                                     vec2mat(legcoef,ngrid) ;
  
  mcurv = -mreconpp ./ (1 + mreconp.^2).^(3/2) ;
  mcurv = mcurv * 337.5 ;
        %  put on scale of diopters

  if isurf == 1 ;
    resdat = mrecon ;
  elseif isurf == 3 ;
    resdat = mcurv ;
  end ;

  rmresid = rmresid - resdat ;


  subplot(3,4,ip) ;
    plot(rgrid,resdat) ;
      title([imagstr ', F=' num2str(mind(ip-1,1)) ...
                     ', L=' num2str(mind(ip-1,2))]) ;
      axis([0 r0 bottom top]) ;
        ssstr = ['% of SS,T = ' num2str(100 * ssoa / sst)] ;
        texty = bottom + .92 * (top - bottom) ;
      text(textx,texty,ssstr) ;  
        ssstr = ['rel. sd = ' num2str(spread)] ;
        texty = bottom + .05 * (top - bottom) ;
      text(textx,texty,ssstr) ;  

end ;

%  ssr is residual sum of squares

disp('    Making Residual plot') ;
subplot(3,4,12) ;
  plot(rgrid,rmresid) ;
    title([imagstr ', res. - ' centstr]) ;
    axis([0 r0 bottom top]) ;
      ssstr = ['% of SS,T = ' num2str(100 * ssr / sst)] ;
      texty = bottom + .9 * (top - bottom) ;
    text(textx,texty,ssstr) ;  



  orient landscape ;
eval(['print -dps \matlab\steve\ps\cd3' num2str(igroup)...
                   num2str(isurf) num2str(idisp) '.ps']) ;





