disp('MATLAB Script File:  arvo1b.m') ;
%
%   For the Cornea Study, this makes ascii files for ARVO poster
%         of April 1997
%   This is for creating ascii files of Fourier 0 mean profiles
%
format compact ;

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



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



for ir0 = 1:2 ;

  if ir0 == 1 ;
    r0 = 4.2 ;    %  Usual radius of analysis
  elseif ir0 == 2 ;
    r0 = 1.7 ;    %  Usual radius of analysis
  end ;


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




  %  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 ;
    if ir0 == 1 ;
      eval(['load \matlab\steve\cornea\fm\fm42' fids '.mat']) ;
    elseif ir0 == 2 ;
      eval(['load \matlab\steve\cornea\fm\fm17' fids '.mat']) ;
    end ;
            %  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') ;



  for ipart = 1:2 ;

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

      if ipart == 1 ;      %  pre-op - others
        if irow == 1 ;     %  Pre Op.
          mfvd = mfv0 ;
        elseif irow == 2 ;     %  1 mo.
          mfvd = mfv1 ;
        elseif irow == 3 ;     %  12 mo.
          mfvd = mfv12 ;
        end ;

      elseif ipart == 2 ;   %  differences
        if irow == 1 ;         %  12 mo. - Pre
          flagp = (vpatind < 400) ;
          mfvd = mfv(:,flagp .* vflag12) - mfv(:,flagp .* vflag0) ;
        elseif irow == 2 ;         %  1 mo. - Pre
          flagp = (vpatind < 500) + (vpatind >= 600) ;
          mfvd = mfv(:,flagp .* vflag1) - mfv(:,flagp .* vflag0) ;
        elseif irow == 3 ;     %  12 mo. - 1 mo.
          flagp = (vpatind < 400) ;
          mfvd = mfv(:,flagp .* vflag12) - mfv(:,flagp .* vflag1) ;
        end ;

      end ;

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

    end ;



    %  WRITE ANSWERS INTO A FILE
    if ir0 == 1 ;
      if ipart == 1 ;
        fstr = '/matlab/steve/arvo42pa.out' ;
        titstr = ['A 4.2 Mean']
      elseif ipart == 2 ;
        fstr = '/matlab/steve/arvo42pb.out' ;
        titstr = ['B 4.2 Diff']
      end ;
    elseif ir0 == 2 ;
      if ipart == 1 ;
        fstr = '/matlab/steve/arvo17pa.out' ;
        titstr = ['A 1.7 Mean']
      elseif ipart == 2 ;
        fstr = '/matlab/steve/arvo17pb.out' ;
        titstr = ['B 1.7 Diff']
      end ;
    end ;
    fid = fopen(fstr,'wt') ;
    
      forstr = '%10.6f,   %10.6f,   %10.6f,   %10.6f\n' ;
              %  Embedded text appears as written (including parens)

    cntbytes = fprintf(fid,forstr,[prgrid mmoutm]') ;

    fclose('all') ;



    
    disp('Making Mean Plot') ;
    subplot(2,2,2*(ipart-1)+ir0) ;
      plot(prgrid,mmoutm) ;
        title(titstr) ;

    
  end ;

end ;

