disp('MATLAB Script File:  cocd1.m') ;
%
%   For the Cornea Study, this does a "curves as data" approach
%       to understanding the variability in the corneal groups.
%   Copied from cosd1.m, which it modifies by changing from an 
%       analysis of the whole surface to just the curves that are
%       the radial components (for a given Fourier order).
%
format compact ;

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

igroup = 2 ;  %  1 - 1 Month
              %  2 - Kerataconic

ifour = 0 ;   %  Fourier Order:   0,1,...,9

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)


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.' ;
elseif igroup == 2 ;
  namestr = '\matlab\steve\cornea\fm\fnameker.asc' ;
  imagstr = 'Ker.' ;
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 ;


  if ifour == 0 ;    %  Then need to get signs right
    vamp = fm(1,:)' ;
  else ;             %  Just use conventional amplitude
    vamp = mamp(ifour+1,:)' ;
          %  Vector of amplitudes for chosen Fourier
  end ;

  mfa = [mfa vamp] ;
          %  put matrix into a column vector and add to right
end ;

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


disp('Working on eigenector decomposition') ;
d = size(mfa,1) ;
n = size(mfa,2) ;
sst = sum(sum(mfa.^2)) ;
%  get mean feature vector, and resids from that
vmean = mean(mfa') ;
vsd = std(mfa') ;
ssm = n * sum(vmean.^2) ;
ssmr = sst - ssm ;
mresid = mfa - vec2mat(vmean',n) ;

%  Get cov's and eigen values
mcov = cov(mresid') ;
          %  Get covariance matrix, transpose, since want "coordinates
          %                      as variables"

[meigvec,meigval] = eig(mcov) ;
veigval = (diag(meigval)) ;
          %  Get eigenvalues and eigenvectors

[veigval,vind] = sort(veigval) ;
meigvec = meigvec(:,vind) ;
          %  Sort to avoid possible ordering problems
  
ev1 = meigvec(:,d) ;
          %  eigenvector with largest eigenvalue
pc1 = mresid' * ev1 ;
          %  vector of projections of data onto ev1
mpc1 = ev1 * pc1' ;
          %  matrix for display
ss1 = sum(sum(mpc1.^2)) ;
ss1r = ssmr - ss1 ;
mresid1 = (mresid - ev1 * pc1') ;

ev2 = meigvec(:,d-1) ;
          %  eigenvector with largest eigenvalue
pc2 = mresid1' * ev2 ;
          %  vector of projections of data onto ev2
mpc2 = ev2 * pc2' ;
          %  matrix for display
ss2 = sum(sum(mpc2.^2)) ;
ss2r = ss1r - ss2 ;
mresid2 = (mresid1 - ev2 * pc2') ;

ev3 = meigvec(:,d-2) ;
          %  eigenvector with largest eigenvalue
pc3 = mresid2' * ev3 ;
          %  vector of projections of data onto ev3
mpc3 = ev3 * pc3' ;
          %  matrix for display
ss3 = sum(sum(mpc3.^2)) ;
ss3r = ss2r - ss3 ;
mresid3 = (mresid2 - ev3 * pc3') ;

ev4 = meigvec(:,d-3) ;
          %  eigenvector with largest eigenvalue
pc4 = mresid3' * ev4 ;
          %  vector of projections of data onto ev4
mpc4 = ev4 * pc4' ;
          %  matrix for display
ss4 = sum(sum(mpc4.^2)) ;
ss4r = ss3r - ss4 ;
mresid4 = (mresid3 - ev4 * pc4') ;


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 Curves') ;
rdat = [] ;
for i = 1:n ;
  legcoef = mfa(:,i) ;

  %  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 ;
    rdat = [rdat mrecon] ;
  elseif isurf == 3 ;
    rdat = [rdat mcurv] ;
  end ;
end ;


disp('Reconstructing mean and sd') ;
legcoef = vmean' ;
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 ;
  rmean = mrecon ;
elseif isurf == 3 ;
  rmean = mcurv ;
end ;
legcoef = vsd' ;
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 ;
  rsd = mrecon ;
elseif isurf == 3 ;
  rsd = mcurv ;
end ;

rmresid = rdat - vec2mat(rmean,n) ;

disp('Reconstructing pc1') ;
rmpc1 = [] ;
for i = 1:n ;
  legcoef = mpc1(:,i) ;
  %  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 ;
    rmpc1 = [rmpc1 mrecon] ;
  elseif isurf == 3 ;
    rmpc1 = [rmpc1 mcurv] ;
  end ;
end ;

rmresid1 = rmresid - rmpc1 ;

disp('Reconstructing pc2') ;
rmpc2 = [] ;
for i = 1:n ;
  legcoef = mpc2(:,i) ;
  %  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 ;
    rmpc2 = [rmpc2 mrecon] ;
  elseif isurf == 3 ;
    rmpc2 = [rmpc2 mcurv] ;
  end ;
end ;

rmresid2 = rmresid1 - rmpc2 ;

disp('Reconstructing pc2') ;
rmpc3 = [] ;
for i = 1:n ;
  legcoef = mpc3(:,i) ;
  %  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 ;
    rmpc3 = [rmpc3 mrecon] ;
  elseif isurf == 3 ;
    rmpc3 = [rmpc3 mcurv] ;
  end ;
end ;

rmresid3 = rmresid2 - rmpc3 ;

disp('Reconstructing pc4') ;
rmpc4 = [] ;
for i = 1:n ;
  legcoef = mpc4(:,i) ;
  %  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 ;
    rmpc4 = [rmpc4 mrecon] ;
  elseif isurf == 3 ;
    rmpc4 = [rmpc4 mcurv] ;
  end ;
end ;

rmresid4 = rmresid3 - rmpc4 ;




%  Make plots
textx = 0.2 ;
if isurf == 1 ;
  bottom1 = -1.2 ;
  top1 = 0 ;
  bottom2 = -.15 ;
  top2 = .15 ;
elseif isurf == 3 ;
  bottom1 = 0 ;
  top1 = 100 ;
  bottom2 = -20 ;
  top2 = 20 ;
end ;
subplot(3,4,1) ;
  plot(rgrid,rdat) ;
    title([imagstr ', Raw Data']) ;
    axis([0 r0 bottom1 top1]) ;
      ssstr = ['SS,T = ' num2str(sst)] ;
      texty = bottom1 + .9 * (top1 - bottom1) ;
    text(textx,texty,ssstr) ;  

  
subplot(3,4,3) ;
  plot(rgrid,rmean,'-',...
       rgrid,rmean + 2 * rsd,':',...
       rgrid,rmean - 2 * rsd,':') ;
    title([imagstr ', Mean +- 2 sd']) ;
    axis([0 r0 bottom1 top1]) ;
      ssstr = ['% of SS,T = ' num2str(100 * ssm / sst)] ;
      texty = bottom1 + .9 * (top1 - bottom1) ;
    text(textx,texty,ssstr) ;  

subplot(3,4,4) ;
  plot(rgrid,rmresid) ;
    title([imagstr ', Res. from Mean']) ;
    axis([0 r0 bottom2 top2]) ;
      ssstr = ['% of SS,T = ' num2str(100 * ssmr / sst)] ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  
  
subplot(3,4,5) ;
  plot(rgrid,rmpc1) ;
    title([imagstr ', PC1']) ;
    axis([0 r0 bottom2 top2]) ;
      ssstr = ['% of SS,MR = ' num2str(100 * ss1 / ssmr)] ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  
    
subplot(3,4,9) ;
  plot(rgrid,rmresid1) ;
    title([imagstr ', Res. from PC1']) ;
    axis([0 r0 bottom2 top2]) ;
      ssstr = ['% of SS,MR = ' num2str(100 * ss1r / ssmr)] ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  

  
subplot(3,4,6) ;
  plot(rgrid,rmpc2) ;
    title([imagstr ', PC2']) ;
    axis([0 r0 bottom2 top2]) ;
      ssstr = ['% of SS,1R = ' num2str(100 * ss2 / ss1r)] ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  
    
subplot(3,4,10) ;
  plot(rgrid,rmresid2) ;
    title([imagstr ', Res. from PC2']) ;
    axis([0 r0 bottom2 top2]) ;
      ssstr = ['% of SS,1R = ' num2str(100 * ss2r / ss1r)] ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  

  
subplot(3,4,7) ;
  plot(rgrid,rmpc3) ;
    title([imagstr ', PC3']) ;
    axis([0 r0 bottom2 top2]) ;
      ssstr = ['% of SS,2R = ' num2str(100 * ss3 / ss2r)] ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  
    
subplot(3,4,11) ;
  plot(rgrid,rmresid3) ;
    title([imagstr ', Res. from PC3']) ;
    axis([0 r0 bottom2 top2]) ;
      ssstr = ['% of SS,2R = ' num2str(100 * ss3r / ss2r)] ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  
  
  
subplot(3,4,8) ;
  plot(rgrid,rmpc4) ;
    title([imagstr ', PC4']) ;
    axis([0 r0 bottom2 top2]) ;
      ssstr = ['% of SS,3R = ' num2str(100 * ss4 / ss3r)] ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  
    
subplot(3,4,12) ;
  plot(rgrid,rmresid4) ;
    title([imagstr ', Res. from PC4']) ;
    axis([0 r0 bottom2 top2]) ;
      ssstr = ['% of SS,3R = ' num2str(100 * ss4r / ss3r)] ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  

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





