disp('MATLAB Script File:  copr4.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 copr3, but now differences are studied
%
format compact ;

autop = 1 ;   %  1 to create a postscript file

ipart = 2 ;   %  1 - study differences pre - 1, pre - 6, pre - 12
              %  2 - study differences 1 - 6, 6 - 12, 1 - 12,

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



r0 = 1.7 ;    %  Smaller 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 ;      %  pre-op - others
  vaxis1 = [0 r0 -20 20] ;
  vaxis2 = [0 r0 -15 15] ;
  vaxis3 = [0 r0 3 7] ;
elseif ipart == 2 ;   %  after surgery differences
  vaxis1 = [0 r0 -20 20] ;
  vaxis2 = [0 r0 -15 15] ;
  vaxis3 = [0 r0 -3 3] ;
end ;

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

  if ipart == 1 ;      %  pre-op - others
    if irow == 1 ;         %  pre-op - 1 mo.
      imagstri = 'PRK, pre - 1 ' ;
      flagp = (vpatind < 500) + (vpatind >= 600) ;
      mfvd = mfv(:,flagp .* vflag0) - mfv(:,flagp .* vflag1) ;
    elseif irow == 2 ;     %  pre-op - 6 mo.
      imagstri = 'PRK, pre - 6 ' ;
      flagp = (vpatind < 200) + ((vpatind >= 400) .* (vpatind < 600)) ;
      mfvd = mfv(:,flagp .* vflag0) - mfv(:,flagp .* vflag6) ;
    elseif irow == 3 ;     %  pre-op - 12 mo.
      imagstri = 'PRK, pre - 12 ' ;
      flagp = (vpatind < 400) ;
      mfvd = mfv(:,flagp .* vflag0) - mfv(:,flagp .* vflag12) ;
    end ;
  elseif ipart == 2 ;   %  after surgery differences
    if irow == 1 ;         %  1 mo. - 6 mo.
      imagstri = 'PRK, 1 - 6 ' ;
      flagp = (vpatind < 200) + ((vpatind >= 400) .* (vpatind < 500)) ;
      mfvd = mfv(:,flagp .* vflag1) - mfv(:,flagp .* vflag6) ;
    elseif irow == 2 ;     %  6 mo. - 12 mo.
      imagstri = 'PRK, 6 - 12 ' ;
      flagp = (vpatind < 200) ;
      mfvd = mfv(:,flagp .* vflag6) - mfv(:,flagp .* vflag12) ;
    elseif irow == 3 ;     %  1 mo. - 12 mo.
      imagstri = 'PRK, 1 - 12 ' ;
      flagp = (vpatind < 400) ;
      mfvd = mfv(:,flagp .* vflag1) - mfv(:,flagp .* vflag12) ;
    end ;
  end ;


  disp('Inverting feature vectors to get profile curves') ;
  moutr = [] ;
  for ipat = 1:size(mfvd,2) ;    %  loop through patients and reconstruct
                                 %        raw profiles
    fmraw = reshape(mfvd(:,ipat),19,7) ;
    [mhr,mhpr,mhppr] = corhpp(fmraw',prgrid,ptgrid,r0) ;
    voutr = -mhppr ./ (1 + mhpr.^2).^(3/2) ;
    moutr = [moutr, voutr * 337.5] ;
    
  end ;


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



  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 ;
    eval(['print \matlab\steve\ps\copr4' num2str(ipart) 'a.ps -dps']) ;
end ;


disp('Any key to continue') ;
pause ;


disp('Making Mean Plot') ;
  subplot(1,1,1) ;
  plot(prgrid,mmoutm) ;
    vh = get(gca,'Children') ;
    set(vh(1),'LineStyle','-.') ;
    set(vh(2),'LineStyle','--') ;
    set(vh(3),'LineStyle','-') ;
    title(['PRK ' zstr ' Means']) ;
    axis(vaxis3) ;
    if ipart == 1 ;
      legend('PO - 1','PO - 6','PO - 12') ;
    elseif ipart == 2 ;
      legend('1 - 6','6 - 12','1 - 12') ;
    end ;

if autop == 1 ;
  orient landscape ;
  eval(['print \matlab\steve\ps\copr4' num2str(ipart) 'b.ps -dps']) ;
end ;




