function mout = coazr(featv) ;
% COAZR Reconstruction of a curvature image from Zernike feature vectors
%      Corneal Topography Analysis System
%   Creates a matrix of curvatures, from a Zernike decomposition
%      Intended to be called by coazrec.m
% Input:
%       featv - Zernike feature vector (in standard order)
% Output:
%       mout - matrix of curvatures, for creating an image

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



%  Set local parameters
%
ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array



%  Set up global variables
%
%  First some general globals
global R0 NMAX IBASIS IRECENT INDATPATH HICURVE ;
%  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
%  IBASIS - 1 for full (triangle), 2 for radial (3 row)
%  IRECENT - To recenter (or not) decompositions
%  INDATPATH - Path with colorfile, cocolor1.dat
%  HICURVE - Gridpoints of Keratron range


%  Calculate grid quantities
%
xgrid = linspace(-R0, R0, ngrid) ;
ygrid = xgrid ;	        %  Vectors of image coordinate values 

[mxgrid, mygrid] = meshgrid(xgrid,ygrid) ;
              		%  Matrices of image coordinate values 
            
mrgrid = sqrt((mxgrid).^2 + (mygrid).^2) ;
              		%  Matrices of image polar coordinate radii
mthetag = atan2(mxgrid,mygrid) ;
              	%  Matrices of image polar coordinate angles (-pi,pi]
vind = find(mthetag < 0) ;
mthetag(vind) = mthetag(vind) + 2 * pi ;
              	%  Matrices of image polar coordinate angles [0,2pi)
          
mcrop = ones(ngrid,ngrid) ;            
  vind = find(mrgrid > R0) ;
              	%  Vector of indices where image is outside circle
mcrop(vind) = nan * ones(length(vind),1) ;
mrgrid(vind) = 0 * ones(length(vind), 1) ;
              	%  Cropping Matrix, multiply by this to cut off images
              	%  outside circle of radius R0.    



ncomap = length(HICURVE) ;
          %  number of rows in the color map ;




%  Construct matrix of curvatures
%
if IBASIS == 1 ;    %  then use full basis
  [mh,mhp,mhpp] = zrecon(featv,NMAX,mrgrid/R0,mthetag,R0) ;
			%  creates vectors of fitted values
			%  for a given data set (ht, deriv, 2deriv)
elseif IBASIS == 2 ;   %  then use radial basis
  [mh,mhp,mhpp] = rrecon(featv,NMAX,mrgrid/R0,mthetag,R0) ;
			%  creates vectors of fitted values
			%  for a given data set (ht, deriv, 2deriv)
end ;


mout = -mhpp ./ (1 + mhp.^2).^(3/2) ;
mout = mout * 337.5 ;
          		%  put on scale of diopters

mout = reshape(mout, [ngrid,ngrid]);
			%  expresses output vector as 'grid matrix'
			%  for plotting purposes


%  Map curvatures to Keratron colors
%
vmout = mout(:) ;
      	   		%  NOTE: interp1 seems to need vector argument
indlo = find(vmout < 0) ;
vmout(indlo) = 0 * ones(length(indlo),1) ;
indup = find(vmout > max(HICURVE)) ;
vmout(indup) = max(HICURVE) * ones(length(indup),1) ;
%  NOTE: interp1 seems to need values within grid

u01 = linspace(0,1,ncomap+1) ;
colcurv = interp1([0; HICURVE], u01, vmout) ;
    	      		%  use this to determine vector of colors representing
          		%  curvature, when plotting.  This transforms curvature 
          		%  to get Keratron colors, when using COMAP as colormap.
          		%  NOTE: interp1 seems to need vector argument
mout = reshape(colcurv,size(mout)) ;  

mout = mout .* mcrop ;
mout = mout' ;




