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


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



%  set local variables
%
ndatin = 2 + 256 * 27 ;




%  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
    %
    RISTR = [INDATPATH infilestr RFILESUF] ;
    HISTR = [INDATPATH infilestr HFILESUF] ;
    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('!!!   Skipping this decomposition') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

    else ;


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

      rsize = length(rdata) ;
      hsize = length(hdata) ;

      if (rsize ~= ndatin) | (hsize ~= ndatin) ;
                             %  then a file is not the right size, so report

        disp('   ') ;
        disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
        disp('!!!   Can''t do Zernike decomposition, ') ;
        disp(['!!!      since a file length is not ' num2str(ndatin) ':']) ;
  
        if (rsize ~= ndatin) ;
          disp(['!!!          Input radius file has length: ' rsize]) ;
        end ;
 
        if (hsize ~= ndatin) ;
          disp(['!!!          Input height file has length: ' hsize]) ;
        end ;

        disp('!!!') ;
        disp('!!!   Skipping this decomposition') ;
        disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

      else ;


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

        rdata = reshape(rdata,256,27) ;
        hdata = reshape(hdata,256,27) ;

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

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



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

        %  first set values outside R0 to missing
        % 
        vor0ind = (rdata > R0) ;
           %  ones where rdata is outside R0
        nor0ind = sum(sum(vor0ind)) ;
        rdata(vor0ind) = zeros(nor0ind,1) ;
        adata(vor0ind) = zeros(nor0ind,1) ;
        hdata(vor0ind) = zeros(nor0ind,1) ;

        %  construct vector of data for "outer boundary"
        %
        rdvmvind = coamdf(rdata) ;
          %  vector of indices before first missing values for rdata
        hdvmvind = coamdf(hdata) ;
          %  vector of indices before first missing values for hdata
        vmvind = min([rdvmvind'; hdvmvind'])' ;
          %  vector of indices before first missing values (overall)
        nonmissr = [] ;
        trimrdata = [] ;
        trimhdata = [] ;
        trimadata = [] ;
        for iang = 1:256 ;
          nonmissr = [nonmissr; rdata(iang,vmvind(iang,1))] ;
          %  radius of outer boundary, for later plotting
          trimrdata = [trimrdata; rdata(iang,1:vmvind(iang,1))'] ;
          trimhdata = [trimhdata; hdata(iang,1:vmvind(iang,1))'] ;
          trimadata = [trimadata; adata(iang,1:vmvind(iang,1))'] ;
          %  trimmed version of data, arranged now as a column
        end ;
          %  take values where are inside R0, 
          %  and are before first non-missing
        rdata = trimrdata ;
        hdata = trimhdata ;
        adata = trimadata ;



        %  Put in form for calling zernike.m
        %  
        rdata = [zeros(256,1); rdata] ;
        adata = [adata(1:256); 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 nonmissr R0 NMAX CRAD CANG']) ;


      end ;    

    end ;

  end ;

end ;










