function coazdec
% COAZDEC, Zernike decomposition
%      Corneal Topography Analysis System, ver. 0.1
%   Creates a file with Zernike Decomposition
%      based on data in the file in the global variable INFILESTR,
%      assumed to reside in the directory
%              in the global variable INDATPATH,
%   in the current Figure Window
%      Intended to be called by coaddm.m
% Output:
%     Matlab data file, with Zernike decomposition in vector form
%            in the directory in the global variable FEATMPATH,


%    Copyright (c) J. S. Marron, N. Locantore 1998
%    Department of Statistics, University of North Carolina



%  Set local parameters
%
icut = 1 ;    %  1 to cut off first two elements in file read
iasciiout = 0 ;     % 1 to write out an ASCII file


%  Set up global variables
%
%  First some general globals
global R0 NMAX CRAD CANG ;
%  R0 - radius of analysis
          %    this gets reloaded when a feature vec is reloaded from file
%  NMAX - maximum number of Zernike columns
          %    this gets reloaded when a feature vec is reloaded from file
%  CRAD - radius of recentering (read in as editable text)
%  CANG - angle of recentering (read in as editable text)

%  Next some globals for file names, etc.
global INFILESTR RFILESUF HFILESUF RISTR HISTR;
%  INFILESTR - Prefix of filename for input data
%  RFILESUF  - radii data file suffix 
%  HFILESUF  - height data file suffix
%  RISTR     - string to evaluate to get radius file	
%  HISTR     - string to evaluate to get height file

%  Next the paths for file control
global INDATPATH FEATMPATH ;
%  INDATPATH - path for input Keratron files
%  FEATMPATH - path for storing Zernike feature matrices



%  Read in data
%
fidr = fopen(RISTR) ;
fidh = fopen(HISTR) ;
          %  assigns a file handle, (read only assumed)
if (fidr == -1) | (fidh == -1) ;      %  then a file not found, so report

  disp('   ') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Can''t do Zernike decomposition,  can''t find file(s):') ;
  
  if (fidr == -1) ;
    disp(['!!!       ' RISTR]) ;
  end ;
 
  if (fidh == -1) ;
    disp(['!!!       ' HISTR]) ;
  end ;

  disp('!!!') ;
  disp('!!!   Terminating raw data display') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

else ;


  %  Continue with data read
  %
  rdata = fscanf(fidr,'%g') ;
            %  formatted reads, assumes single column
  hdata = fscanf(fidh,'%g') ;
  fclose(fidr) ;
  fclose(fidh) ;


  %  Put in matrix form, and do preliminary adjustments
  %
  if icut == 1 ;     %  then need to cut off first two elements
    rdata = rdata(3:length(rdata)) ;
    hdata = hdata(3:length(hdata)) ;
  end ;

  adatasingle = linspace(0, (1 - 1/256) * 2 * pi, 256)' ;
  adata = kron(ones(27,1), adatasingle) ;
          %  Corresponding matrix of angles.

  rdata = rdata(257:length(rdata)) ;
  adata = adata(257:length(adata)) ;
  hdata = hdata(257:length(hdata)) ;
          %  cutoff first columns of all zeros.


  %  remove inappropriate values
  %
  %  first cutoff missings in radii files
  vind = find(rdata == 0) ;
           %  vector of indices where have a missing
  rdata(vind) = [] ;

  adata(vind) = [] ;
  hdata(vind) = [] ;
           %  Deletes that part of the vectors

  %  now cutoff missings in height files
  vind = find(hdata == 0) ;
           %  vector of indices where have a missing
  rdata(vind) = [] ;
  adata(vind) = [] ;
  hdata(vind) = [] ;
           %  Deletes that part of the vectors
  %disp(['After removing missings, have ' num2str(length(rdata)) ' values']) ;

  %  now cutoff the values outside R0
  vind = find(rdata > R0) ;
           %  vector of indices where outside
  rdata(vind) = [] ;
  adata(vind) = [] ;
  hdata(vind) = [] ;
           %  Deletes that part of the vectors
  %disp(['After removing pts > R0, have ' num2str(length(rdata)) ' values']) ;


  %  Put in form for calling zernike.m
  %  
  rdata = [zeros(256,1); rdata] ;
  adata = [adatasingle; adata] ;
  hdata = [zeros(256,1); hdata] ;

  rdata = rdata / R0 ;
          %  normalize to Zernike radius 1


  %  Do least square fitting
  %
  designm = zernike(NMAX,R0,rdata,adata) ;
          %  construct Zernike design matrix 
  betahat = (designm' * designm) \ (designm' * hdata) ;
          %  do ordinary least squares fit


  outfilestr = [FEATMPATH INFILESTR '.mat'] ;  
  eval(['save ' outfilestr ' betahat R0 NMAX CRAD CANG']) ;
  
  if iasciiout == 1 ;    %  then also create an ASCII file

      outfilestr = [FEATMPATH INFILESTR '.asc'] ;  
    fid = fopen(outfilestr,'wt') ;
          %  'wt' is for "delete contents of this file and open 
          %               for writing" (with 't' for "text").

      forstr = '%12.10\n' ;
    cntbytes = fprintf(fid,forstr,betahat) ;

  end ;

end ;







