disp('MATLAB Script File:  cocd4.m') ;
%
%   For the Cornea Study, this does a "curves as data" approach
%       to understanding the variability in the corneal groups.
%   Copied from cocd4.m, which it modifies by using full
%       groups, by doing a "full decomposition", and by allowing
%       looking at variability, either "raw", or "about the mean",
%       or about a "central image".
%
format compact ;

idel = 1 ;    %  1 to delete outliers (currently only for groups 7,9)

autop = 1 ;   %  1 to create a postscript file

ipause = 1 ;  %  0 to just make all pics, etc., with no pauses
              %  1 to pause between each plot

igroup = 9 ;  %  1 - Normal       (Central Image = "All K")
              %  2 - Ker Suspect  (Central Image = "All K")
              %  3 - Kerataconic  (Central Image = "All K")
              %  4 - Pre          (Central Image = "All P")
              %  5 - 1 Mo.        (Central Image = "All P")
              %  6 - 1 Yr.        (Central Image = "All P")
              %  7 - Pre - 1 Mo.  (Central Image not inter.)
              %  8 - Pre - 1 Yr.  (Central Image not inter.)
              %  9 - 1 Mo - 1 Yr  (Central Image not inter.)

ictype = 0 ;   %  0 - no centering
               %  1 - center about "central image" for this group
               %  2 - center about mean for this group

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

izoom = 2 ;   %  zoom in (has effect for height maps only)
              %  0 - full scale showing all heights
              %  1 - first zoom
              %  2 - higher magnification zoom
              %          (very useful for difference curves)

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


if igroup == 1 ;
  namestr = '\matlab\steve\cornea\fm\fnnor.asc' ;
  imagstr = 'Norm' ;
  icent = 119 ;    %  cofmap1.m index for Centering Feature Matrix
              %  119 - Mean All K
  centstr = '- all K' ;
  nincol = 1 ;
elseif igroup == 2 ;
  namestr = '\matlab\steve\cornea\fm\fnksus.asc' ;
  imagstr = 'Ksus' ;
  icent = 119 ;    %  cofmap1.m index for Centering Feature Matrix
              %  119 - Mean All K
  centstr = '- all K' ;
  nincol = 1 ;
elseif igroup == 3 ;
  namestr = '\matlab\steve\cornea\fm\fnker.asc' ;
  imagstr = 'Ker' ;
  icent = 119 ;    %  cofmap1.m index for Centering Feature Matrix
              %  119 - Mean All K
  centstr = '- all K' ;
  nincol = 1 ;
elseif igroup == 4 ;
  namestr = '\matlab\steve\cornea\fm\fnpre.asc' ;
  imagstr = 'Pre' ;
  icent = 122 ;    %  cofmap1.m index for Centering Feature Matrix
              %  122 - Mean all P
  centstr = '- all P' ;
  nincol = 1 ;
elseif igroup == 5 ;
  namestr = '\matlab\steve\cornea\fm\fn1mo.asc' ;
  imagstr = '1 Mo.' ;
  icent = 122 ;    %  cofmap1.m index for Centering Feature Matrix
              %  122 - Mean all P
  centstr = '- all P' ;
  nincol = 1 ;
elseif igroup == 6 ;
  namestr = '\matlab\steve\cornea\fm\fn1yr.asc' ;
  imagstr = '1 Yr.' ;
  icent = 122 ;    %  cofmap1.m index for Centering Feature Matrix
              %  122 - Mean all P
  centstr = '- all P' ;
  nincol = 1 ;
elseif igroup == 7 ;
  namestr = '\matlab\steve\cornea\fm\fnpre1mo.asc' ;
  imagstr = 'Pre - 1Mo' ;
  icent = 0 ;   %  no need to center differences
  centstr = '' ;
  nincol = 2 ;
elseif igroup == 8 ;
  namestr = '\matlab\steve\cornea\fm\fnpre1yr.asc' ;
  imagstr = 'Pre - 1Yr' ;
  icent = 0 ;   %  no need to center differences
  centstr = '' ;
  nincol = 2 ;
elseif igroup == 9 ;
  namestr = '\matlab\steve\cornea\fm\fn1mo1yr.asc' ;
  imagstr = '1Mo - 1Yr' ;
  icent = 0 ;   %  no need to center differences
  centstr = '' ;
  nincol = 2 ;
end ;


r0 = 4.2 ;    %  Usual radius of analysis

ngrid = 100 ; 

nr = 10 ;
nc = 7 ;



disp('Loading feature matrices for this group') ;
fid = fopen(namestr,'r') ;
if nincol == 1 ;
  mfindx = fscanf(fid,'%4d') 
elseif nincol == 2 ;
  mfindx = fscanf(fid,'%4d %4d') ;
  mfindx = reshape(mfindx,2,length(mfindx)/2)' ;
end ;

%  Used this for checking
%mfindx
%pause ;


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

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

  if nincol == 2 ;    %  Then need to subtract second feature matrix
    fm1 = fm ;
    if mfindx(i,2) < 10 ;
      eval(['load ' prestr '000' num2str(mfindx(i,2))]) ;
    elseif mfindx(i,2) < 100 ;
      eval(['load ' prestr '00' num2str(mfindx(i,2))]) ;
    elseif mfindx(i,2) < 1000 ;
      eval(['load ' prestr '0' num2str(mfindx(i,2))]) ;
    else ;
      eval(['load ' prestr num2str(mfindx(i,2))]) ;
    end ;
          %  This creates a number of variables, especially fm,
          %  a 19 x 7 "feature matrix"
    fm = fm1 - fm ;
  end ;

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


if igroup == 7 | igroup == 9 ;
  if idel == 1 ;
    if igroup == 7 ;
      mfa = [mfa(:,1:26) mfa(:,28:42)] ;
    elseif igroup == 9 ;
      mfa = [mfa(:,1:26) mfa(:,28)] ;
    end ;
  end ;
end ;


if  ictype > 0  &  nincol == 1  ;
  disp('Do Centering') ;
  if ictype == 1 ;   %  Then do "central image" centering
    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(mfindx,1)) ;
  elseif ictype == 2 ;  %  Then Center at Mean

    mfa = mfa - vec2mat(mean(mfa')',size(mfindx,1)) ;

    centstr = '- mean' ;
  end ;
else ;  %  No centering is done
  centstr = '' ;
end ;



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

disp(['joint 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('   Working on 1st plot') ;
figure(1) ;
mout = zeros(ngrid,n) ;
vssf = zeros(10,1) ;
for ifour = 1:10 ;    %  Loop through Fourier levels

  eval(['mout' num2str(ifour) ' = zeros(ngrid,n) ;']) ;
  for i = 1:n ;    %  Loop through individuals
    legcoef = mfa((0:10:60)+ifour,i) ;
    vssf(ifour) = vssf(ifour) + sum(legcoef .^ 2) ;

    %  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 ;
      mout(:,i) = mout(:,i) + mrecon ;
      eval(['mout' num2str(ifour) '(:,i) = mrecon ;']) ;
    elseif isurf == 3 ;
      mout(:,i) = mout(:,i) + mcurv ;
      eval(['mout' num2str(ifour) '(:,i) = mcurv ;']) ;
    end ;
 
  end ;
end ;



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

%  Reset Graphics to get proper blanks
clf ;

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

subplot(3,4,12) ;
  bar(0:9,100 * vssf / sst) ;
    title([imagstr ', Raw ' centstr]) ;
    axis([-.5 9.5 0 100]) ;
    text(-.2,92,'% of SS,T by F') ;

for ifour = 1:10 ;    %  Loop through Fourier levels
  subplot(3,4,ifour+1) ;
    eval(['plot(rgrid,mout' num2str(ifour) ') ;']) ;
      title([imagstr ', F' num2str(ifour-1) ' ' centstr]) ;
      axis([0 r0 bottom top]) ;
        ssstr = ['% of SS,T = ' num2str(100 * vssf(ifour) / sst)] ;
        texty = bottom + .92 * (top - bottom) ;
      text(textx,texty,ssstr) ;  
end ;

if autop == 1 ;
    orient landscape ;
  eval(['print -dps \matlab\steve\ps\cd4' num2str(igroup)...
               num2str(isurf) num2str(ictype) 'a.ps']) ;
end ;


if ipause == 1 ;
  disp('Hit any key to do decompositions of Fourier pieces') ;
  pause ;
end ;

mssfl = [] ;
for ifour = 1:10 ;    %  Loop through Fourier levels
  disp(['    Working on Fourier ' num2str(ifour - 1)]) ;
  figure(2) ;

  mout = zeros(ngrid,n) ;
  vssfl = zeros(1,7) ;
  for ileg = 1:7 ;    %  Loop through Legendre orders
    eval(['mout' num2str(ileg) ' = zeros(ngrid,n) ;']) ;
  end ;

  for i = 1:n ;    %  Loop through individuals
    legcoef = mfa((0:10:60)+ifour,i) ;
    vssfl = vssfl + (legcoef .^ 2)' ;

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

    if isurf == 1 ;
      mout(:,i) = sum(mrecon)' ;
    elseif isurf == 3 ;
      mout(:,i) = sum(mcurv)' ;
    end ;
 
    for ileg = 1:7 ;    %  Loop through Legendre orders
      if isurf == 1 ;
        eval(['mout' num2str(ileg) '(:,i) = mrecon(ileg,:)'' ;']) ;
      elseif isurf == 3 ;
        eval(['mout' num2str(ileg) '(:,i) = mcurv(ileg,:)'' ;']) ;
      end ;
    end ;

  end ;



  subplot(3,3,1) ;
    plot(rgrid,mout) ;
      title([imagstr ', F' num2str(ifour-1) ' ' centstr]) ;
      axis([0 r0 bottom top]) ;
        ssstr = ['% of SS,T = ' num2str(100 * vssf(ifour) / sst)] ;
        texty = bottom + .9 * (top - bottom) ;
      text(textx,texty,ssstr) ;  

  subplot(3,3,9) ;
    bar(0:6,100 * vssfl / vssf(ifour)) ;
      title([imagstr ', F' num2str(ifour-1) ' ' centstr]) ;
      axis([-.5 6.5 0 100]) ;
      text(-.2,92,'% of SS,F by L') ;


  for ileg = 1:7 ;    %  Loop through Legendre orders
    subplot(3,3,ileg+1) ;
      eval(['plot(rgrid,mout' num2str(ileg) ') ;']) ;
       title([imagstr ', F' num2str(ifour-1) ' L' num2str(ileg-1) ...
                                             centstr]) ;
        axis([0 r0 bottom top]) ;
          ssstr = ['% of SS,T = ' num2str(100 * vssfl(ileg) / sst)] ;
          texty = bottom + .92 * (top - bottom) ;
        text(textx,texty,ssstr) ;  
          ssstr = ['% of SS,F = ' ...
                            num2str(100 * vssfl(ileg) / vssf(ifour))] ;
          texty = bottom + .04 * (top - bottom) ;
        text(textx,texty,ssstr) ;  

  end ;

  if autop == 1 ;
      orient landscape ;
    eval(['print -dps \matlab\steve\ps\cd4' num2str(igroup)...
                 num2str(isurf) num2str(ictype) num2str(ifour-1) '.ps']) ;
  end ;


  if ipause == 1 ;
    disp('        Any key for next one') ;
    pause ;
  end ;

  mssfl = [mssfl; vssfl] ;

end ;


%  Used these for Checking
%disp('Check these are 0') ;
%abs(sum(sum(mssfl)) - sst) 
%sum(abs(vssf - sum(mssfl')')) 
%pause ;

vssl = sum(mssfl) ;


%  Write results to a file
  foutstr = ['\matlab\steve\ps\cd4' num2str(igroup) ...
                       num2str(ictype) '.out'] ;
fid = fopen(foutstr,'wt') ;
    titlstr1 = 'Output from the MATLAB Script cocd4.m,   ' ;
    titlstr1 = [titlstr1,date] ;
    titlstr2 = '  Does decomposition of Sums of Squares for Cornea Groups' ;
  cntbytes = fprintf(fid,'%1s\n',titlstr1) ;
  cntbytes = fprintf(fid,'%1s\n\n\n',titlstr2) ;

  hdrstr1 = ['For the Group: ' imagstr ',    Centered by: ' centstr] ;
  cntbytes = fprintf(fid,'%1s\n\n\n',hdrstr1) ;

  hdrstr2 = ['                         % of Total Sum of Squares'] ;
  hdrstr3 = ['                L0      L1      L2      L3      L4'] ;
  hdrstr3 = [hdrstr3, '      L5      L6'] ;
  cntbytes = fprintf(fid,'%1s\n',hdrstr2) ;
  cntbytes = fprintf(fid,'%1s\n',hdrstr3) ;

    fstr =   '       all    %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n\n' ;
    vout = 100 * vssl / sst ;
  cntbytes = fprintf(fid,fstr,vout) ;
  for iout = 0:9 ;
      fstr = ['F' num2str(iout)] ;
      fstr = [fstr, '  %7.3f   %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n'] ;
      vout = 100 * [vssf(iout+1) mssfl(iout+1,:)] / sst ;
    cntbytes = fprintf(fid,fstr,vout) ;
  end ;


  hdrstr2 = ['                         % of Fourier Sums of Squares'] ;
  hdrstr3 = ['                L0      L1      L2      L3      L4'] ;
  hdrstr3 = [hdrstr3, '      L5      L6'] ;
  cntbytes = fprintf(fid,'\n\n\n%1s\n',hdrstr2) ;
  cntbytes = fprintf(fid,'%1s\n',hdrstr3) ;

    fstr =   '       all    %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n\n' ;
    vout = 100 * vssl / sst ;
  cntbytes = fprintf(fid,fstr,vout) ;
  for iout = 0:9 ;
      fstr = ['F' num2str(iout)] ;
      fstr = [fstr, '  %7.3f   %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n'] ;
      vout = 100 * [vssf(iout+1) mssfl(iout+1,:)] / vssf(iout+1) ;
    cntbytes = fprintf(fid,fstr,vout) ;
  end ;


  hdrstr2 = ['                         % of Legendre Sums of Squares'] ;
  hdrstr3 = ['                L0      L1      L2      L3      L4'] ;
  hdrstr3 = [hdrstr3, '      L5      L6'] ;
  cntbytes = fprintf(fid,'\n\n\n%1s\n',hdrstr2) ;
  cntbytes = fprintf(fid,'%1s\n',hdrstr3) ;

    fstr =   '       all    %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n\n' ;
    vout = 100 * (vssl ./ vssl) ;
  cntbytes = fprintf(fid,fstr,vout) ;
  for iout = 0:9 ;
      fstr = ['F' num2str(iout)] ;
      fstr = [fstr, '  %7.3f   %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n'] ;
      vout = 100 * [(vssf(iout+1)/sst) (mssfl(iout+1,:) ./ vssl)] ;
    cntbytes = fprintf(fid,fstr,vout) ;
  end ;


fclose('all') ;


