disp('MATLAB Script File: coleg11.m') ;
%
%   For the Cornea Study, this does twiddles with Legendre's
%       including interpretability
%   This looks at cumulative construction of Legendre's
%       2nd Generation, using complete data sets
%
%
format compact ;

ipart = 1 ;   %  1 - do Legendre's for first 4 Fouriers

autop = 1 ;   %  1 - for auto print

igroup = 21    %  1 - Repeatability Study (all)
              %  20 - PRK Study (Pre Op)
              %  21 - PRK Study (1 mo)
              %  26 - PRK Study (6 mo)
              %  212 - PRK Study (12 mo)
              %  300 - Kerataconic (normal)
              %  301 - Kerataconic (susp map)
              %  302 - Kerataconic (unilat)
              %  303 - Kerataconic (clinical)

ivpat = 3 ;

if ivpat == 1 ;
  vpat = [100; 104; 107; 402] ;    % interesting for igroup = 21 
elseif ivpat == 2 ;
  vpat = [411; 102; 103; 204] ;    % expect good fit, got inverse u
elseif ivpat == 3 ;
  vpat = [205; 216; 217; 209] ;    % expect poor fit, got u
end ;




r0 = 4.2 ;    %  Usual radius of analysis
vaxis1 = [0 r0 25 65] ;

ngrid = 27 ;

rgrid = linspace(0,r0,ngrid)' ;


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




%  stuff for loading feature matrices
if igroup == 1 ;  %  Work with Repeatability population
  imagstr = 'Rep''t ' ;
  disp(['Loading feature matrices for ' imagstr]) ;
  %  First load in file info from file rptkey.dat, 
  %  which was modified from Jinting's MATLAB script normind.m
  %  (everything but the data was edited out).
  %    1st col is patient id
  %    2nd col is file id (add last 4 digits to "rpfm" to get file suffix)
  %    3rd col is "eye"
  %    4th col is pupil distance
  %    5th col is pupil direction
  fid = fopen('\matlab\steve\cornea\fm\rptkey.dat','r') ;
    frmtstr = '%g %g %g %g %g' ;
      mfilein = fscanf(fid,frmtstr,[5,inf]) ;
      mfilein = mfilein' ;
            %  since data rows are read in as columns
    fclose(fid) ;

  nimag = size(mfilein,1) ;
  vpatind = mfilein(:,1) ;
            %  vector of patient id's
  vfileid = mfilein(:,2) ;
            %  vector of file id's
  veyeind = mfilein(:,3) ;
            %  vector of eye indices
  mfilein = 0 ;
            %  won't be using this again


  %  Loop through feature matrix files, and combine into big data matrix
  mfv = [] ;
  for iimag = 1:nimag ;
    eval(['load \matlab\steve\cornea\fm\rpfm' ...
                             num2str(vfileid(iimag)) '.mat']) ;
            %  This load Jinting's previously saved stuff
            %  which is a number of variables, especially fm,
            %  a 19 x 7 "feature matrix"
    mfv = [mfv fm(:)] ;
            %  put matrix into a column vector and add to right
            %  Including the feature matrix for this image
  end ;


elseif igroup < 300 ;   %  work with prk population
  disp('Working on PRK Population ') ;
  %  First load in file info from file prkkey.dat, 
  %  which was modified from fcfilesn.txt
  %  (everything but the data was edited out).
  %    1st col is patient id
  %    2nd col is Time (0,1,6,12 months)
  %    3nd col is file id (add last 4 digits to "fm42" to get file suffix)
  %    4th col is "eye"
  %    5th col is pupil distance
  %    6th col is pupil direction
  fid = fopen('\matlab\steve\cornea\fm\prkkey.dat','r') ;
    frmtstr = '%g %g %g %g %g %g' ;
      mfilein = fscanf(fid,frmtstr,[6,inf]) ;
      mfilein = mfilein' ;
            %  since data rows are read in as columns
    fclose(fid) ;

  vpatind = mfilein(:,1) ;
            %  vector of patient id's
  vtime = mfilein(:,2) ;
            %  vector of time's
  vfileid = mfilein(:,3) ;
            %  vector of file id's
  veyeind = mfilein(:,4) ;
            %  vector of eye indices
  mfilein = 0 ;
            %  won't be using this again


  if igroup == 20 ;
    imagstr = 'PRK Pre-op' ;
    disp(['Loading feature matrices for ' imagstr]) ;
    vflag = (vtime == 0) ;
  elseif igroup == 21 ;
    imagstr = 'PRK 1 Mo' ;
    disp(['Loading feature matrices for ' imagstr]) ;
    vflag = (vtime == 1) ;
  elseif igroup == 26 ;
    imagstr = 'PRK 6 Mo' ;
    disp(['Loading feature matrices for ' imagstr]) ;
    vflag = (vtime == 6) ;
  elseif igroup == 212 ;
    imagstr = 'PRK 12 Mo' ;
    disp(['Loading feature matrices for ' imagstr]) ;
    vflag = (vtime == 12) ;
  end ;

  vpatind = vpatind(vflag) ;
  vfileid = vfileid(vflag) ;
  veyeind = veyeind(vflag) ;
  nimag = length(vpatind) ;


  %  Loop through feature matrix files, and combine into big data matrix
  mfv = [] ;
  for iimag = 1:nimag ;
    fid = vfileid(iimag) ;
    if fid < 100 ;
      fids = ['00' num2str(fid)] ;
    elseif fid < 1000 ;
      fids = ['0' num2str(fid)] ;
    else ;
      fids = num2str(fid) ;
    end ;
    eval(['load \matlab\steve\cornea\fm\fm42' fids '.mat']) ;
            %  This load Jinting's previously saved stuff
            %  which is a number of variables, especially fm,
            %  a 19 x 7 "feature matrix"
    mfv = [mfv fm(:)] ;
            %  put matrix into a column vector and add to right
            %  Including the feature matrix for this image
  end ;


else ;    %  work with Keratoconus group
  disp('Working on Ker Population ') ;
  %  First load in file info from file kerkey.dat, 
  %  which was modified from fcfilesn.txt
  %  (everything but the data was edited out).
  %    1st col is patient id
  %    2nd col is Ker type Code
  %    3nd col is file id (add last 4 digits to "fm42" to get file suffix)
  %    4th col is "eye"
  %    5th col is pupil distance
  %    6th col is pupil direction
  fid = fopen('\matlab\steve\cornea\fm\kerkey.dat','r') ;
    frmtstr = '%g %g %g %g %g %g' ;
      mfilein = fscanf(fid,frmtstr,[6,inf]) ;
      mfilein = mfilein' ;
            %  since data rows are read in as columns
    fclose(fid) ;

  vpatind = mfilein(:,1) ;
            %  vector of patient id's
  vkercode = mfilein(:,2) ;
            %  vector of time's
  vfileid = mfilein(:,3) ;
            %  vector of file id's
  veyeind = mfilein(:,4) ;
            %  vector of eye indices
  mfilein = 0 ;
            %  won't be using this again


  if igroup == 300 ;
    imagstr = 'Ker Normal' ;
    disp(['Loading feature matrices for ' imagstr]) ;
    vflag = (vkercode == 0) ;
  elseif igroup == 301 ;
    imagstr = 'Ker sus map' ;
    disp(['Loading feature matrices for ' imagstr]) ;
    vflag = (vkercode == 1) ;
  elseif igroup == 302 ;
    imagstr = 'Ker unilat' ;
    disp(['Loading feature matrices for ' imagstr]) ;
    vflag = (vkercode == 2) ;
  elseif igroup == 303 ;
    imagstr = 'Ker clinical' ;
    disp(['Loading feature matrices for ' imagstr]) ;
    vflag = (vkercode == 3) ;
  end ;

  vpatind = vpatind(vflag) ;
  vfileid = vfileid(vflag) ;
  veyeind = veyeind(vflag) ;
  nimag = length(vpatind) ;


  %  Loop through feature matrix files, and combine into big data matrix
  mfv = [] ;
  for iimag = 1:nimag ;
    fid = vfileid(iimag) ;
    if fid < 100 ;
      fids = ['00' num2str(fid)] ;
    elseif fid < 1000 ;
      fids = ['0' num2str(fid)] ;
    else ;
      fids = num2str(fid) ;
    end ;
    eval(['load \matlab\steve\cornea\fm\fm42' fids '.mat']) ;
            %  This load Jinting's previously saved stuff
            %  which is a number of variables, especially fm,
            %  a 19 x 7 "feature matrix"
    mfv = [mfv fm(:)] ;
            %  put matrix into a column vector and add to right
            %  Including the feature matrix for this image
  end ;


end ;

fclose('all') ;


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




sphht = spherad - sqrt(spherad.^2 - rgrid.^2) ;
%plot(rgrid,sphht) ;



if ipart == 1 ;    %  make plots using 0th Fourier coefficients

  for ipat = 1:4 ;

    flag = (vpatind == vpat(ipat)) ;
    if sum(flag) ~= 1 ;
      disp('!!!   Careful, wrong number of patients   !!!') ;
    end ;
    fv = mfv(:,flag) ;
    mfeat = reshape(fv,19,7) ;
    mfeat = mfeat' ;

    legcoef = mfeat(:,1) ;
        %  this is Fourier 0

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


    subplot(2,2,ipat) ;
      plot(rgrid,mcurv) ;
        title(['Cum L''s for F0 ' num2str(vpat(ipat)) ', Curv']) ;
          legend('l=0','l=1','l=2','l=3','l=4','l=5','l=6',-1) ;
        axis(vaxis1) ;
  end ;

  if autop ~= 0 ;
    eval(['print -dps le11' num2str(igroup) num2str(ivpat) '.ps']) ;
  end ;

elseif ipart == 2 ;


end ;

