disp('Running MATLAB script file vire1.m') ;
%

%    Reads already stored Zernike Decomps, and saves as
%    ASCII files, for e.g. SAS analysis
%    Also does first analyses
%



ipart = 2 ;      %  1 - Create Ascii files
                 %  2 - Some Sum of Square Analyses

ipat = 0 ;       %  0 - loop through all
                 %  1 - Hickson
                 %  2 - Meadows
                 %  3 - Tesoro 
                 %  4 - Yi     

inewfig = 0 ;    %  0 - all pics in one figure window
                 %  1 - make a new figure window for each pic

iR0 = 2 ;        %  1 - R0 = 4.2
                 %  2 - R0 = 3.0


%  Initialize parameters
%
if ipat == 0 ;     %  then loop through all patients
  ips = 1 ;
  ipe = 4 ;
else ;             %  then just do input patients
  ips = ipat ;
  ipe = ipat ;
end ;

ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array
xgrid = linspace(-R0, R0, ngrid) ;
ygrid = xgrid ;	        %  Vectors of image coordinate values 

global R0 NMAX IBASIS IRECENT INDATPATH HICURVE ;
NMAX = 10 ;
IBASIS = 1 ;
IRECENT = 0 ;

if iR0 == 1 ;
  outsuf = 'a' ;
  R0 = 4.2 ;
elseif iR0 == 2 ;
  outsuf = 'b' ;  
  R0 = 3.0 ;
end ;

INDATPATH = '\matlab52\cornean\indata\' ;
          %  path for input Keratron files
[COMAP, HICURVE] = coacmap ;





if ipart == 1 ;    %  then create big ASCII file with results


  %  Do appropriate file reads
  %
  mout = [] ;
  for ip = ips:ipe ;

    if ip == 1 ;    
      patname = 'Hickson' ;
      featmpath = '\matlab52\cornean\zfm\virehick\' ;
      listr = [featmpath 'hickall.lst'] ;
    elseif ip == 2 ;    
      patname = 'Meadows' ;
      featmpath = '\matlab52\cornean\zfm\viremead\' ;
      listr = [featmpath 'meadall.lst'] ;
    elseif ip == 3 ;    
      patname = 'Tesoro' ;
      featmpath = '\matlab52\cornean\zfm\vireteso\' ;
      listr = [featmpath 'tesoall.lst'] ;
    elseif ip == 4 ;    
      patname = 'Yi' ;
      featmpath = '\matlab52\cornean\zfm\vireyi\' ;
      listr = [featmpath 'yiall.lst'] ;
    end ;



    %  read files, using lines from   coalwzr.m
    %
    fid = fopen(listr) ;
            %  assigns a file handle, (read only assumed)

    if (fid == -1) ;      %  then file not found, so report
      disp('   ') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      disp('!!!   Warning from vire1.m:   can''t find file:') ;
      disp(['!!!       ' listr]) ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    end ;


    %  Read in list of filenames
    %
    minfilestr = fscanf(fid,'%s',[8,inf]) ;
              %  formatted read, assumes single column
    fclose(fid) ;

    minfilestr = minfilestr' ;
            %  each row corresponds to one filename


    %  loop through files (one for each image)
    %
    ie = size(minfilestr,1) ;
    for i = 1:ie ;

      infilestr = minfilestr(i,:) ;

      disp(['  Loading ' num2str(i) ' of ' num2str(ie) ...
                               ' for ' patname ', file: ' infilestr ]) ;


      zmistr = [featmpath infilestr '.mat'] ;
            %  full path and filename 
      eval(['load ' zmistr]) ;
            %  loads this Zernike feature matrix and parameter values

      vout = betahat' ;
            %  turn column vector of coefficients into a row vector

      filenum = str2num(infilestr(1:6)) ;
      filenum = filenum + (str2num(infilestr(8))) / 10 ;
            %  get first 6 numbers as an integer, 
            %  then add last digit after decimal

      vout = [filenum ip vout] ;
            %  insert ASCII numerical file number, and patient number


      mout = [mout; vout] ;
            %  tack on new row to main output matrix


    end ;

  end ;



  %  Save as ASCII file, using lines from gptable.m
  %
  fid = fopen(['\matlab52\cornean\vire1' outsuf '.out'],'wt') ;
          %  'wt' is for "delete contents of this file and open 
          %               for writing" (with 't' for "text").
          %  Use 'ay' to append to what is already there, 'rt' for reading.



  forstr = '%6.1f %1.0f' ;
  for iformat = 1:length(betahat) ;
    forstr = [forstr ' %8.6f'] ;
  end ;
  forstr = [forstr '\n'] ;

  cntbytes = fprintf(fid,forstr,mout') ;
          %  Careful here, because this reads down columns first,
          %  So feed it the transpose of the usual form
  fclose(fid) ;


  disp(' ') ;
  disp(' ') ;
  disp(['Finished constructing ASCII file: vire1' outsuf '.out']) ;
  disp(['Now need to edit in Day, Time, R-L, ']) ;
  disp(['and save as: vire1' outsuf '.in']) ;



elseif ipart == 2 ;    % then do some decompositions of sums of squares


  %  First read in data file, editted from file created above
  %
  fid = fopen(['\matlab52\cornean\vire1' outsuf '.in'],'rt') ;
          %  'rt' is for "read only" and "text"

  forstr = '%f %f %f %f %f' ;
          %  format for reading:
          %      file no.
          %      pat. no.
          %      day no.
          %      session no.
          %      eye no.
  ninpar = 5 ;
  ncoeff = 66 ;
  for iformat = 1:ncoeff ;
    forstr = [forstr ' %f'] ;
  end ;

  [inmat,cnt] = fscanf(fid,forstr) ;
  fclose(fid) ;

  ncol = ninpar + ncoeff ;
  nrow = cnt / ncol ;

  inmat = reshape(inmat,ncol,nrow)' ;
          %   since reshape works column wise

%  inmat(1:10,:)
%  inmat(nrow,:)
          %  used for testing read in



  grandmean = mean(inmat(:,(ninpar+1):ncol)) ;
          %  vector of coefficents of mean of everything


  %  Main loop throught patients
  %
  for ip = ips:ipe ;

    if ip == 1 ;    
      patname = 'Hick' ;
    elseif ip == 2 ;    
      patname = 'Mead' ;
    elseif ip == 3 ;    
      patname = 'Teso' ;
    elseif ip == 4 ;    
      patname = 'Yi' ;
    end ;

    pflag = (inmat(:,2) == ip) ;
          %  2nd column is patient number
    pdat = inmat(pflag,:) ;
          %  matrix of data for this patient
    pmean = mean(pdat(:,(ninpar+1):ncol)) ;
          %  vector of coefficents of mean for this patient


    for ieye = 1:2 ;      %  loop through eyes

      peflag = (pdat(:,5) == ieye) ;
          %  5th column is eye number
      pedat = pdat(peflag,:) ;
          %  matrix of data for this patient, this eye
      pemean = mean(pedat(:,(ninpar+1):ncol)) ;
          %  vector of coefficents of mean for this patient, this eye


      if inewfig == 0 ;
        figure(1) ;
        clf ;
      elseif inewfig == 1 ;
        figure ;
      end ;


      for iday = 1:3 ;      %  loop though days

        pedflag = (pedat(:,3) == iday) ;
          %  3rd column is day number
        peddat = pedat(pedflag,:) ;
          %  matrix of data for this patient, this eye, this day
        pedmean = mean(peddat(:,(ninpar+1):ncol)) ;
          %  vector of coefficents of mean for patient, eye, day


        for isess = 1:3 ;     %  loop through sessions

          pedsflag = (peddat(:,4) == isess) ;
            %  4th column is session number

          nimag = sum(pedsflag) ;
          dispstr = ['Av''ging ' num2str(nimag) ' images, for ' patname ...
                        ', eye ' num2str(ieye) ', day ' num2str(iday) ...
                        ', session ' num2str(isess)] ;
          disp(dispstr) ;

          if nimag > 0 ;
            pedsdat = peddat(pedsflag,:) ;
              %  matrix of data for this patient, this eye, this day, this sess
            pedsmean = mean(pedsdat(:,(ninpar+1):ncol)) ;
              %  vector of coefficents of mean for patient, eye, day, sess

            subplot(4,6,6*(iday-1)+isess+1) ;
            titstr = [patname ' E_' num2str(ieye) ' D_' num2str(iday) ...
                              ' S_' num2str(isess) ' Avg.'] ;

            %  Create image with Keratron Color scheme, and put in current axes
            %
            betahat = pedsmean' ;
            imagout = coazr(betahat) ;
            pcolor(xgrid,ygrid,imagout) ;
              shading flat ;
              axis('square') ;
              title(titstr) ;
                caxis([0 1]) ;
              colormap(COMAP) ;

            hold on ;
              tx = -.95 * R0 ;
              ty = .7 * R0 ;
              iss = nimag * sum((pedsmean - grandmean).^2) ;
              tss = sum(sum((pedsdat(:,(ninpar+1):ncol) - ...
                                       vec2mat(grandmean,nimag)).^2)) ;
              ssstr = num2str(100 * iss / tss) ;
              tstr = [ssstr '% of GMR'] ;
              text(tx,ty,tstr) ;
            hold off ;


          end;

        end ;


        %  Make marginal plots for Days
        %
        subplot(4,6,6*(iday-1)+5) ;
         titstr = [patname ' E_' num2str(ieye) ' D_' num2str(iday) ...
                           ' Avg.'] ;

         %  Create image with Keratron Color scheme, and put in current axes
         %
         betahat = pedmean' ;
         imagout = coazr(betahat) ;
         pcolor(xgrid,ygrid,imagout) ;
           shading flat ;
           axis('square') ;
           title(titstr) ;
             caxis([0 1]) ;
           colormap(COMAP) ;

         hold on ;
           tx = -.95 * R0 ;
           ty = .7 * R0 ;
           nimag = size(peddat,1) ;
           iss = nimag * sum((pedmean - grandmean).^2) ;
           tss = sum(sum((peddat(:,(ninpar+1):ncol) - ...
                                 vec2mat(grandmean,nimag)).^2)) ;
           ssstr = num2str(100 * iss / tss) ;
           tstr = [ssstr '% of GMR'] ;
           text(tx,ty,tstr) ;
         hold off ;


      end ;



      %  Make marginal plots for Sessions
      %
      for isess = 1:3 ;     %  loop through sessions

        pesflag = (pedat(:,4) == isess) ;
            %  4th column is session number

        pesdat = pedat(pesflag,:) ;
              %  matrix of data for this patient, this eye, this sess
        pesmean = mean(pesdat(:,(ninpar+1):ncol)) ;
              %  vector of coefficents of mean for patient, eye, sess

        subplot(4,6,6*3+isess+1) ;
        titstr = [patname ' E_' num2str(ieye) ...
                              ' S_' num2str(isess) ' Avg.'] ;

        %  Create image with Keratron Color scheme, and put in current axes
        %
        betahat = pesmean' ;
        imagout = coazr(betahat) ;
        pcolor(xgrid,ygrid,imagout) ;
          shading flat ;
          axis('square') ;
          title(titstr) ;
            caxis([0 1]) ;
          colormap(COMAP) ;

        hold on ;
          tx = -.95 * R0 ;
          ty = .7 * R0 ;
          nimag = size(pesdat,1) ;
          iss = nimag * sum((pesmean - grandmean).^2) ;
          tss = sum(sum((pesdat(:,(ninpar+1):ncol) - ...
                                   vec2mat(grandmean,nimag)).^2)) ;
          ssstr = num2str(100 * iss / tss) ;
          tstr = [ssstr '% of GMR'] ;
          text(tx,ty,tstr) ;
        hold off ;

      end ;


      %  Make overall Mean plot
      %
      subplot(4,6,6*3+5) ;
      titstr = [patname ' E_' num2str(ieye) ' Avg.'] ;

      %  Create image with Keratron Color scheme, and put in current axes
      %
      betahat = pemean' ;
      imagout = coazr(betahat) ;
      pcolor(xgrid,ygrid,imagout) ;
        shading flat ;
        axis('square') ;
        title(titstr) ;
          caxis([0 1]) ;
        colormap(COMAP) ;

      hold on ;
        tx = -.95 * R0 ;
        ty = .7 * R0 ;
        nimag = size(pedat,1) ;
        iss = nimag * sum((pemean - grandmean).^2) ;
        tss = sum(sum((pedat(:,(ninpar+1):ncol) - ...
                                 vec2mat(grandmean,nimag)).^2)) ;
        ssstr = num2str(100 * iss / tss) ;
        tstr = [ssstr '% of GMR'] ;
        text(tx,ty,tstr) ;
      hold off ;


        orient landscape ;
      eval(['print -dpsc \matlab52\cornean\vr1' outsuf num2str(ip) ...
                             num2str(ieye) '.ps']) ;


    end ;



  end ;


end ;


