disp('MATLAB Script File:  copr2.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 takes a "transformation" approach, where the curves
%       are "registered" via linear transformation
%
format compact ;

autop = 1 ;   %  1 to create a postscript file

ipart = 1 ;   %  1 - for each of 1,6,12 mo, show raw profiles,
              %               registered profiles and both means
              %  2 - Overlay means
              %  3 - Find cases with immediate descent from center

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

itrans = 1 ;  %  1 - just scales x axis


r0 = 4.2 ;    %  Usual radius of analysis

nrgrid = 43 ;
prgrid = linspace(0,4.2,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 "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
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\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 ;

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 4.2 25 65] ;
  vaxis2 = [0 4.2 35 51] ;

  moutm = [] ;
  mouttm = [] ;
  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 ;
    n = size(mfvi,2) ;
          %  number of patients
    d = length(prgrid) ;
          %  number of grid points

    disp('Inverting feature vectors to get profile curves') ;
    moutr = [] ;
    for ipat = 1:n ;    %  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) ;
    voutm = -mhppr ./ (1 + mhpr.^2).^(3/2) ;
    voutm = voutm * 337.5 ;
    
    disp('Making Raw Plot') ;
    subplot(3,3,3*(irow-1)+1) ;
      plot(prgrid,moutr) ;
        title([imagstri zstr ' Raw']) ;
        axis(vaxis1) ;

    
    if itrans == 1 ;
      transtr = ' x-scal' ;

      [vmaxht,vmaxloc] = max(moutr) ;
      vmaxloc = prgrid(vmaxloc) ;
      maxlmean = mean(vmaxloc) ;
           %  mean of maximum locations
      vmaxloc = vmaxloc' ;

      vmag = maxlmean ./ vmaxloc ;
           %  vector of magnification factors to make means same
      mtrgrid = vec2mat(prgrid,n) .*  vec2mat(vmag,d) ;
           %  matrix of transformed rgrids

      %  Loop through vectors, to form transformed mean
      vouttm = zeros(d,1) ;
      tmdenom = zeros(d,1) ;
      for ipat = 1:n ;
        flag = prgrid > mtrgrid(d,ipat) ;
        if sum(flag) == 0 ;    %  then mag data extends beyond end
          ilast = d ;
        else ;         %  then have "blank spots"
          [temp,ilast] = max(flag) ;
          ilast = ilast - 1 ;
        end ;

        vouttm(1:ilast) = vouttm(1:ilast) + ...
               interp1(mtrgrid(:,ipat),moutr(:,ipat),prgrid(1:ilast)) ;
        tmdenom(1:ilast) = tmdenom(1:ilast) + ones(ilast,1) ;
      end ;

      vouttm = vouttm ./ tmdenom ;

    end ;


    disp('Making Tranformed Plot') ;
    subplot(3,3,3*(irow-1)+2) ;
      plot(mtrgrid,moutr) ;
        title([imagstri transtr]) ;
        axis(vaxis1) ;

    
    disp('Making Mean Plots') ;
    subplot(3,3,3*(irow-1)+3) ;
      plot(prgrid,voutm,'--',prgrid,vouttm,'-') ;
        title([imagstri zstr ' Means']) ;
%        legend('raw','x-scaled') ;
        axis(vaxis2) ;
    
    
    moutm = [moutm voutm] ;
    mouttm = [mouttm vouttm] ;


  end ;

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



  disp('Any key for combined mean plots') ;
  pause ;

    disp('Making Big Mean Plots') ;
    subplot(2,1,1) ;
      plot(prgrid,moutm) ;
        vh = get(gca,'Children') ;
        set(vh(1),'LineStyle','-.') ;
        set(vh(2),'LineStyle','--') ;
        set(vh(3),'LineStyle','-') ;
        title(['PRK Reg Study, ' zstr ' Raw Means']) ;
        axis(vaxis2) ;
        legend('1 mo','6 mo','12 mo') ;

    subplot(2,1,2) ;
      plot(prgrid,mouttm) ;
        vh = get(gca,'Children') ;
        set(vh(1),'LineStyle','-.') ;
        set(vh(2),'LineStyle','--') ;
        set(vh(3),'LineStyle','-') ;
        title(['PRK Reg Study, ' zstr ' X-Scaled Means']) ;
        axis(vaxis2) ;
        legend('1 mo','6 mo','12 mo') ;


  if autop == 1 ;
    orient tall ;
      print \matlab\steve\ps\copr21b.ps -dps ;
  end ;



elseif ipart == 2 ;      %  then 


elseif ipart == 3 ;      %  then find cases with decent from center

  diary copr2.out ;
  disp('Patient Numbers, where decreases near 0') ;
  diary off ;

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

    if irow == 1 ;     %  1 mo.
      mfvi = mfv1 ;
      imagstri = imagstr1 ;
      vpatindi = vpatind1 ;
    elseif irow == 2 ;     %  6 mo.
      mfvi = mfv6 ;
      imagstri = imagstr6 ;
      vpatindi = vpatind6 ;
    elseif irow == 3 ;     %  12 mo.
      mfvi = mfv12 ;
      imagstri = imagstr12 ;
      vpatindi = vpatind12 ;
    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 ;
    
    vflagd = (moutr(2,:) - moutr(1,:)) < 0 ;
          %  one where derivative at 0 negative

    disp('Making Decreasing Plot') ;
    subplot(3,3,3*(irow-1)+1) ;
      plot(prgrid,moutr(:,vflagd)) ;
        title([imagstri zstr ' Dec.']) ;
        axis(vaxis1) ;

    disp('Making Increasing Plot') ;
    subplot(3,3,3*(irow-1)+2) ;
      plot(prgrid,moutr(:,1-vflagd)) ;
        title([imagstri zstr ' Inc.']) ;
        axis(vaxis1) ;

   diary on ;
     disp(imagstri) ;
     vpatindi(vflagd) 
  diary off ;

  end ;


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




end ;



