function coalwzd
% COALWZD, List-Wise Zernike Decomposition
%      Corneal Topography Analysis System, ver. 0.1
%   Creates several files with Zernike Decompositions
%      based on data from files in the global variable LISTFILESTR
%      All of the listed files are assumed to reside in the directory
%              in the global variable INDATPATH,
%   in the current Figure Window
%      Intended to be called by coaddm.m
%      Reads list of files from LISTFILESTR
% 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
okread = 0 ;    %  1 to  indicate successful file read
nang = 256; 	% number of semi-meridians read per radial ring
%  set local variables
ndatin = 2 + 256 * 27 ;


%  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)
          %    this gets reloaded when a feature vec is reloaded from file
% CANG, angle of recentering (read in as editable text)
          %    this gets reloaded when a feature vec is reloaded from file



%  Next some globals for file names, etc.
global RFILESUF HFILESUF RISTR HISTR LISTR ;
%  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
%  LISTR     - string for inputting list file

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

global RFILESUF HFILESUF ;
        %  File Suffixes for data files


%  Currently only do this with no recentering, 
%  so set these to 0 for correct save
%
CANG = 0 ;
CRAD = 0 ;


%  Read in filename
%
% LISTR, (read-in file string, set in advance in coapam.m)
fid = fopen(LISTR) ;
          %  assigns a file handle, (read only assumed)
if (fid == -1) ;      %  then file not found, so report

  disp('   ') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Can''t do Zernike Decomps, because can''t find file:') ;
  disp(['!!!       ' LISTR]) ;
  disp('!!!') ;
  disp('!!!   Terminating List Wise Zernike Decomp.') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

else ;


  %  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 filenames
  %
  nfile = size(minfilestr,1) ;
          %  number of files = rows of input matrix
  for ifile = 1:nfile ;

    infilestr = minfilestr(ifile,:) ;

    disp(['        Decomposing data file ' infilestr]) ;

    %  Read in data
    %
    okread = 0 ;	% initialize successful read flag on each pass
    RISTR = [INDATPATH infilestr RFILESUF] ;
    HISTR = [INDATPATH infilestr HFILESUF] ;
    fidr = fopen(RISTR, 'r', 'l') ;
	%  IEEE format for binary files
    fidh = fopen(HISTR, 'r', 'l') ;
          %  assigns a file handle, (read only assumed)
    if (fidr == -1) | (fidh == -1) ;      %  then a file not found, so reformat
      if (upper(RFILESUF) == RFILESUF);
         RFILESUF = lower(RFILESUF);		% handles case sensitivity
         HFILESUF = lower(HFILESUF);		% .xlb vs. .XLB
      else				
         RFILESUF = upper(RFILESUF);
         HFILESUF = upper(HFILESUF);
      end;
      newristr = [INDATPATH infilestr RFILESUF];
      newhistr = [INDATPATH infilestr HFILESUF];
      fidr = fopen(newristr, 'r', 'l');
      fidh = fopen(newhistr, 'r', 'l');
           %   try binary open with different suffix
    end;

    if (fidr == -1) | (fidh == -1) ;    % then a file not found, so try ascii
      newrsuf = upper([RFILESUF(1:3) 'a']);
      newhsuf = upper([HFILESUF(1:3) 'a']);
      newristr = [INDATPATH infilestr newrsuf];
      newhistr = [INDATPATH infilestr newhsuf];
      fidr = fopen(newristr);
      fidh = fopen(newhistr);
	   	 %   new file handles for ascii format

      if (fidr == -1) | (fidh == -1) ;    % then a file not found, so reformat
        newrsuf = lower(newrsuf);		% handles case sensitivity
        newhsuf = lower(newhsuf);		% .xla vs. .XLA
        newristr = [INDATPATH infilestr newrsuf];
        newhistr = [INDATPATH infilestr newhsuf];
        fidr = fopen(newristr);
        fidh = fopen(newhistr);
            %   try ascii open with different suffix
      end;

   	if (fidr == -1) | (fidh == -1) ;   % then file not found, so give error
     	  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('Skipping this decomposition.') ;
     	  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!');

   	else ;

  	%  Continue with ascii data read, do some data preprocessing
  	%
   	  okread = 1; 
    	  disp('Performing ASCII file read for the input data.');
    	  disp('');
    	  rdata = fscanf(fidr,'%g') ;
               		%  formatted reads, assumes single column
     	  hdata = fscanf(fidh,'%g') ;
     	  fclose(fidr) ;
     	  fclose(fidh) ;
     
     	  if rdata(2) == 0 ;		% (no flag on front of data)
	     nring = (length(rdata)/nang) + 1;
	       if nring - floor(nring) > 0
	      	disp('!!! File not read in properly, dimensions are off  !!!');
	      	okread = 0;
          end;
     	  else
          nring = rdata(2) + 1;
          rdata(1:2) = [];
	       hdata(1:2) = [];
          if nring - floor(nring) > 0
	         disp('!!! File not read in properly, dimensions are off  !!!');
            okread = 0;
          end;  
        end;
      end; 

    else ;

   	%  Continue with binary data read, do some data preprocessing
   	%
   	okread = 1;
  	   disp('Performing binary file read for the input data.');
   	disp('');
   	rstatus = fseek(fidr, 2, -1);
   	rdata = fread(fidr, 'float');
   	fclose(fidr);

   	hstatus = fseek(fidh, 2, -1);
   	hdata = fread(fidh, 'float');
   	fclose(fidh);

   	if rdata(1) == 0 ;		% (no flag on front of data)
   	  nring = (length(rdata)/nang) + 1;
     	  if nring - floor(nring) > 0 ; 
    	     disp('!!! File not read in properly, dimensions are off  !!!');	
	        okread = 0;
     	  end;
   	else
     	  nring = rdata(1) + 1;
     	  rdata(1) = [];
     	  hdata(1) = [];
     	  if nring - floor(nring) > 0 ;
	     disp('!!! File not read in properly, dimensions are off  !!!');
	     okread = 0;
     	  end;  
   	end;
    end;

    if okread == 1; 	% indicates successful file reading

      rdata = reshape(rdata,nang,nring) ;
      hdata = reshape(hdata,nang,nring) ;

      adata = linspace(0, (1 - 1/nang) * 2 * pi, nang) ;
      hold1 = adata;
		% holds a vector of angles for the tracing of boundary
      adata = adata' * ones(1,nring) ;
      	        %  Corresponding matrix of angles.

      rdata = rdata(:,2:nring) ;
      adata = adata(:,2:nring) ;
      hdata = hdata(:,2:nring) ;
          %  cutoff first columns of all zeros

        %  Now get boundary for non-missing values, and restrict data
        %

        %  construct vector of data for inner and outer boundary
        %
   	[rinner, rdvmvind] = coamdf(rdata, nring-1, nang) ;
   	       %  vector of indices before first missing values for rdata
   	[hinner, hdvmvind] = coamdf(hdata, nring-1, nang) ;
  	        %  vector of indices before first missing values for hdata
   	vmvind = min([rdvmvind'; hdvmvind'])' ;
  	        %  vector of indices before first missing values (exterior)
   	vinner = max([rinner'; hinner'])' ;
   	       %  vector of indices after first missing values (interior)

   	nonmissr = [] ;
   	inmissr  = [] ;
   	trimrdata = [] ;
   	trimhdata = [] ;
   	trimadata = [] ;
    for iang = 1:nang ;
  	   inmissr  = [inmissr; rdata(iang, vinner(iang,1))] ;
	 	 %  radius of inner boundary, for later plotting
	   nonmissr = [nonmissr; rdata(iang,vmvind(iang,1))] ;
   	       %  radius of outer boundary, for later plotting
    	trimrdata = [trimrdata; rdata(iang,vinner(iang,1):vmvind(iang,1))'] ;
    	trimhdata = [trimhdata; hdata(iang,vinner(iang,1):vmvind(iang,1))'] ;
    	trimadata = [trimadata; adata(iang,vinner(iang,1):vmvind(iang,1))'] ;
     	     %  trimmed version of data, arranged now as one long column
    end ;
  	        %  take values before first non-missing
      rdata = trimrdata ;
      hdata = trimhdata ;
      adata = trimadata ;

    	%  Reattach data from zero radius
      %  
      rdata = [zeros(nang,1); rdata] ;
      adata = [adata(1:nang); adata] ;
      hdata = [zeros(nang,1); hdata] ;

      %  now cutoff the values outside R0
      %
      vind = find(rdata > 1.05*R0) ;
           %  vector of indices where outside
      rdata(vind) = [] ;
      adata(vind) = [] ;
      hdata(vind) = [] ;
           %  Deletes that part of the vectors


   	%  Do least square fitting
   	%

      rdata = rdata / R0 ;
          %  normalize to Zernike radius 1
      designm = zernike(NMAX,R0,rdata,adata) ;
          %  construct Zernike design matrix 
      betahat = (designm' * designm) \ (designm' * hdata) ;
          %  do ordinary least squares fit
	   radsym = 0;
      	  %  tells us that this is a normal decomposition

      outfilestr = [FEATMPATH infilestr '.mat'] ;  
      eval(['save ' outfilestr ' betahat nonmissr ' ... 
	                   'inmissr radsym R0 NMAX CRAD CANG']);
   	disp(['Finished decomposition for file ' infilestr]) ;
	   disp('');
    end;

  end ;    

end ;




