disp('MATLAB Script File:  copr3.m') ;
%
%   For the Cornea Study, 2nd Generation, using complete data sets
%   This does visual analysis of the "regression - progression"
%       issue for the PRK study
%       but only for the 0th Fourier coefficients 
%       (i.e. in spherical direction only)
%     Color maps are now replaced by curve overlays
%   This is similar to pr1, but now the analysis region is 
%       reduced to r_0 = 1.7 (instead of usual 4.2)
%
format compact ;

autop = 1 ;   %  1 to create a postscript file

ipart = 3 ;   %  1 - for each of 1,6,12 mo, show raw profile, mean and resid
              %  2 - overlay 3 means 1,6,12
              %  3 - similar for Pre-op
              %  4 - overlay 4 means 0,1,6,12

izero = 1 ;   %  0 - Use full feature matrix
              %  1 - Use only Fourier 0 terms



r0 = 1.7 ;    %  Usual radius of analysis

nrgrid = 35 ;
prgrid = linspace(0,r0,nrgrid)' ;
          %  grid for plotting profiles          
ptgrid = zeros(nrgrid,1) ;



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 "fm17" 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
nimag = length(vpatind) ;


imagstr0 = 'PRK Pre-op' ;
vflag0 = (vtime == 0) ;
vpatind0 = vpatind(vflag0) ;
vfileid0 = vfileid(vflag0) ;
veyeind0 = veyeind(vflag0) ;
nimag0 = length(vpatind0) ;

imagstr1 = 'PRK 1 Mo' ;
vflag1 = (vtime == 1) ;
vpatind1 = vpatind(vflag1) ;
vfileid1 = vfileid(vflag1) ;
veyeind1 = veyeind(vflag1) ;
nimag1 = length(vpatind1) ;

imagstr6 = 'PRK 6 Mo' ;
vflag6 = (vtime == 6) ;
vpatind6 = vpatind(vflag6) ;
vfileid6 = vfileid(vflag6) ;
veyeind6 = veyeind(vflag6) ;
nimag6 = length(vpatind6) ;

imagstr12 = 'PRK 12 Mo' ;
vflag12 = (vtime == 12) ;
vpatind12 = vpatind(vflag12) ;
vfileid12 = vfileid(vflag12) ;
veyeind12 = veyeind(vflag12) ;
nimag12 = length(vpatind12) ;




%  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\fm17' 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 ;

fclose('all') ;




%  Now zero out terms in the feature matrix,
if izero == 0 ;     %  Use Full Feature Matrix
  zstr = ' Full FM, ' ;
elseif izero == 1 ;     %  Use only Fourier 0 terms
  zstr = ' Only F0, ' ;
  mask = [ones(1,7); zeros(18,7)] ; 
          %  1's where want to keep coefficients in 19x7 form
  mask = vec2mat(mask(:),size(mfv,2)) ;
          %  First put into a long column, then extend horizontally
  mfv = mfv .* mask ;
end ;



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


mfv0 = mfv(:,vflag0) ;
mfv1 = mfv(:,vflag1) ;
mfv6 = mfv(:,vflag6) ;
mfv12 = mfv(:,vflag12) ;

vmean = mean(mfv') ;
vmean0 = mean(mfv0') ;
vmean1 = mean(mfv1') ;
vmean6 = mean(mfv6') ;
vmean12 = mean(mfv12') ;



if ipart == 1 ;      %  then show 3 means and residuals

  vaxis1 = [0 r0 25 50] ;
  vaxis2 = [0 r0 -10 10] ;


  for irow = 1:3 ;    %  loop through 3 rows of output plot: 1,6,12

    if irow == 1 ;     %  1 mo.
      mfvi = mfv1 ;
      vmeani = vmean1 ;
      imagstri = imagstr1 ;
    elseif irow == 2 ;     %  6 mo.
      mfvi = mfv6 ;
      vmeani = vmean6 ;
      imagstri = imagstr6 ;
    elseif irow == 3 ;     %  12 mo.
      mfvi = mfv12 ;
      vmeani = vmean12 ;
      imagstri = imagstr12 ;
    end ;

    disp('Inverting feature vectors to get profile curves') ;
    moutr = [] ;
    for ipat = 1:size(mfvi,2) ;    %  loop through patients and reconstruct
                                  %        raw profiles
      fmraw = reshape(mfvi(:,ipat),19,7) ;
      [mhr,mhpr,mhppr] = corhpp(fmraw',prgrid,ptgrid,r0) ;
      voutr = -mhppr ./ (1 + mhpr.^2).^(3/2) ;
      moutr = [moutr, voutr * 337.5] ;
      
    end ;
    
    fmmean = reshape(vmeani,19,7) ;
    [mhr,mhpr,mhppr] = corhpp(fmmean',prgrid,ptgrid,r0) ;
    moutm = -mhppr ./ (1 + mhpr.^2).^(3/2) ;
    moutm = moutm * 337.5 ;
    
    disp('Making Raw Plot') ;
    subplot(3,3,3*(irow-1)+1) ;
      plot(prgrid,moutr) ;
        title([imagstri zstr ' Raw']) ;
        axis(vaxis1) ;
    
    
    disp('Making Mean Plot') ;
    subplot(3,3,3*(irow-1)+2) ;
      plot(prgrid,moutm) ;
        title([imagstri zstr ' Mean']) ;
        axis(vaxis1) ;
    
    disp('Making Mean Resid. Plot') ;
    resid = moutr - vec2mat(moutm,size(moutr,2)) ;
    subplot(3,3,3*(irow-1)+3) ;
      plot(prgrid,resid) ;
        title([imagstri zstr ' Mean R']) ;
        axis(vaxis2) ;
    

  end ;


  if autop == 1 ;
    orient landscape ;
      print \matlab\steve\ps\copr31.ps -dps ;
  end ;



elseif ipart == 2 ;      %  then overlay 3 means

  vaxis1 = [0 r0 38 42] ;

  subplot(1,1,1) ;

  mmoutm = [] ;
  for irow = 1:3 ;    %  loop through 3 rows of output plot: 1,6,12

    if irow == 1 ;     %  1 mo.
      mfvi = mfv1 ;
      vmeani = vmean1 ;
      imagstri = imagstr1 ;
    elseif irow == 2 ;     %  6 mo.
      mfvi = mfv6 ;
      vmeani = vmean6 ;
      imagstri = imagstr6 ;
    elseif irow == 3 ;     %  12 mo.
      mfvi = mfv12 ;
      vmeani = vmean12 ;
      imagstri = imagstr12 ;
    end ;

    fmmean = reshape(vmeani,19,7) ;
    [mhr,mhpr,mhppr] = corhpp(fmmean',prgrid,ptgrid,r0) ;
    moutm = -mhppr ./ (1 + mhpr.^2).^(3/2) ;
    mmoutm = [mmoutm moutm * 337.5] ;

  end ;


  disp('Making Mean Plot') ;
    plot(prgrid,mmoutm) ;
      vh = get(gca,'Children') ;
      set(vh(1),'LineStyle','-.') ;
      set(vh(2),'LineStyle','--') ;
      set(vh(3),'LineStyle','-') ;
      title([imagstri zstr ' Mean']) ;
      axis(vaxis1) ;
      legend('1 mo','6 mo','12 mo') ;


  if autop == 1 ;
    orient landscape ;
      print \matlab\steve\ps\copr32.ps -dps ;
  end ;


elseif ipart == 3 ;      %  then show pre-op means and residuals

  vaxis1 = [0 r0 25 50] ;
  vaxis2 = [0 r0 -10 10] ;



  mfvi = mfv0 ;
  vmeani = vmean0 ;
  imagstri = imagstr0 ;

  disp('Inverting feature vectors to get profile curves') ;
  moutr = [] ;
  for ipat = 1:size(mfvi,2) ;    %  loop through patients and reconstruct
                                  %        raw profiles
    fmraw = reshape(mfvi(:,ipat),19,7) ;
    [mhr,mhpr,mhppr] = corhpp(fmraw',prgrid,ptgrid,r0) ;
    voutr = -mhppr ./ (1 + mhpr.^2).^(3/2) ;
    moutr = [moutr, voutr * 337.5] ;
      
  end ;
    
  fmmean = reshape(vmeani,19,7) ;
  [mhr,mhpr,mhppr] = corhpp(fmmean',prgrid,ptgrid,r0) ;
  moutm = -mhppr ./ (1 + mhpr.^2).^(3/2) ;
  moutm = moutm * 337.5 ;
    
  disp('Making Raw Plot') ;
  subplot(3,3,1) ;
    plot(prgrid,moutr) ;
      title([imagstri zstr ' Raw']) ;
      axis(vaxis1) ;
    
    
  disp('Making Mean Plot') ;
  subplot(3,3,2) ;
    plot(prgrid,moutm) ;
      title([imagstri zstr ' Mean']) ;
      axis(vaxis1) ;
  
  disp('Making Mean Resid. Plot') ;
  resid = moutr - vec2mat(moutm,size(moutr,2)) ;
  subplot(3,3,3) ;
    plot(prgrid,resid) ;
      title([imagstri zstr ' Mean R']) ;
      axis(vaxis2) ;
    


  if autop == 1 ;
    orient landscape ;
      print \matlab\steve\ps\copr33.ps -dps ;
  end ;



elseif ipart == 4 ;      %  then overlay 4 means

  vaxis1 = [0 r0 38 45] ;

  subplot(1,1,1) ;

  mmoutm = [] ;
  for irow = 0:3 ;    %  loop through 4 rows of output plot: 0,1,6,12

    if irow == 0 ;     %  Pre op
      mfvi = mfv0 ;
      vmeani = vmean0 ;
      imagstri = imagstr0 ;
    elseif irow == 1 ;     %  1 mo.
      mfvi = mfv1 ;
      vmeani = vmean1 ;
      imagstri = imagstr1 ;
    elseif irow == 2 ;     %  6 mo.
      mfvi = mfv6 ;
      vmeani = vmean6 ;
      imagstri = imagstr6 ;
    elseif irow == 3 ;     %  12 mo.
      mfvi = mfv12 ;
      vmeani = vmean12 ;
      imagstri = imagstr12 ;
    end ;

    fmmean = reshape(vmeani,19,7) ;
    [mhr,mhpr,mhppr] = corhpp(fmmean',prgrid,ptgrid,r0) ;
    moutm = -mhppr ./ (1 + mhpr.^2).^(3/2) ;
    mmoutm = [mmoutm moutm * 337.5] ;

  end ;


  disp('Making Mean Plot') ;
    plot(prgrid,mmoutm) ;
      vh = get(gca,'Children') ;
      set(vh(1),'LineStyle','-.') ;
      set(vh(2),'LineStyle','--') ;
      set(vh(3),'LineStyle','-') ;
      set(vh(4),'LineStyle','*') ;
      title(['PRK ' zstr ' Mean']) ;
      axis(vaxis1) ;
      legend('Pre-Opt','1 mo','6 mo','12 mo') ;


  if autop == 1 ;
    orient landscape ;
      print \matlab\steve\ps\copr34.ps -dps ;
  end ;


end ;



