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
%
okread = 0 ;    %  1 to  indicate successful file read
nang = 256; 	% number of semi-meridians read per radial ring
iasciiout = 1 ;     % 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)
          %    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 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, '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('!!!   Terminating 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 we have original data with all missing values thrown out.
%###  Need to transform this data, and then get rid of all values
%###	which are outside of R0.
%###  also need to transform the boundary region for the plot

%  #########################################

  if (~(CRAD==0))
	[xbd, ybd] = pol2cart(hold1', nonmissr);
			% takes boundary
	angc = CANG * pi/180;		%% put on radian scale

	[xp, yp] = pol2cart(angc, CRAD);
				% converts polar pupil coords to cartesian
	[xdata, ydata] = pol2cart(adata, rdata);
				% converts polar data coords to cartesian
	away = sqrt((xdata - xp).^2 + (ydata - yp).^2);
				% distance from each point to pupil center
	[away2, index] = sort(away);	
			% gives positions of minimum distances
	min1 = [xdata(index(1)), ydata(index(1)), hdata(index(1))];

	j = 2;
	samept = 1;

	while samept ==1;	%% checks that the x/y are not the same (origin)
		min2 = [xdata(index(j)), ydata(index(j)), hdata(index(j))];
		if(min1(1)==min2(1) & min1(2)==min2(2))
			j = j + 1;
		else
			samept = 0;
			v1 = min2 - min1;	% first vec to determine plane
		end
	end

	colinear = 1;		% error check for 3 points in plane  
	j = 3;			% start with third smallest distance

	while colinear == 1
		min3 = [xdata(index(j)), ydata(index(j)), hdata(index(j))];
					% try j-th closest to make plane
		v2 = min3 - min1;	% second vector to determine plane
		if rank([v1;v2]) == 2
			colinear = 0 ;	% found independent vector
		else 
			j = j + 1;	% found colinear vector
		end
	end

	det1 = [1 0 0]; det2 = [0 1 0]; det3= [0 0 1];
	u1 = det([det1; v1; v2]);
	u2 = det([det2; v1; v2]);
	u3 = det([det3; v1; v2]);
	u = [u1; u2; u3];		
					% normal vector (yields ax+by+cz = ?)
	const = min1 * u;		% constant for ? above
	if u3 ~= 0
		u = -u/u3;		% rescale so that Z-direction is +
		const = -const/u3;	%  "
	end

	%% now normalize the vector

	u = u / norm(u);
	const = const / norm(u);

%% find approximate height at pupil center
%% solve for c in the ax+by+cz=d equation found above
%% take care that min(3)'s third value can be zero!

	hold2 = const - min3(1:2)*u(1:2);
	zp = hold2 / min3(3);

	if isnan(zp)
		hold3 = const - min2(1:2)*u(1:2);
		zp = hold3 / min2(3);		% need another check here?
	end;

	pupil = [xp yp zp];		% coordinates of pupil center point.

%% find pupil center axes
%% 

	e1 = [1 0 0]'; 
	e2 = [0 1 0]';

	e1hat = e1 - (u'*e1)*u;
	e1hat = e1hat / norm(e1hat);		% normalized "nasal axis"

	e2hat = e2 - (u'*e2)*u - (e1hat'*e2)*e1hat;
	e2hat = e2hat / norm(e2hat);		% normalized "superior axis"

%% take data and transform the heights H to heights H-twiddle
%%

	xstar = xp - xdata;		
	ystar = yp - ydata;
	zstar = zp - hdata;

	htwiddle = [xstar ystar zstar]*u;	% get transformed heights

%% do the same for the boundary points
%%
	xbdstar = xp - xbd;
	ybdstar = yp - ybd;
	zbdstar = zp - zeros(size(xbd));

%% get transformed x- and y-axes.
%% 

	hatmatrix = eye(3) - u*u';
	ztwiddle = hatmatrix * [-xstar -ystar -zstar]';	
					% transformed positions of data

	xtwiddle = ztwiddle' * e1hat;
	ytwiddle = ztwiddle' * e2hat;
					% x, y values in transformed space

	[adata, rdata] = cart2pol(xtwiddle, ytwiddle);
					% expresses values in polar coords
	hdata = htwiddle;
	disp(['Using recentered data for decomposition']);

%% do same for boundary points
%%
	zbdtwid = hatmatrix * [-xbdstar -ybdstar -zbdstar]';
	xbdtwid = zbdtwid' * e1hat;
	ybdtwid = zbdtwid' * e2hat;

	[hold4, nonmissr] = cart2pol(xbdtwid, ybdtwid);

	end
%  #########################################
%  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
%%
%% set nonmissr to correct value
	nonmissr = min(nonmissr, R0);

  %  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']) ;
  
  if iasciiout == 1 ;    %  then also create an ASCII file

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

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

    fclose(fid) ;

  end ;


end ;







