function coarawd
% COARAWD, Display of Raw Keratron Data Files, in Cornean
%      Corneal Topography Analysis System, ver. 0.1
%   Creates a plot of the raw curvature image,
%      for the file in the global variable INFILESTR,
%      assumed to reside in the directory
%              in the global variable INDATPATH,
%   in the current Figure Window
%   Note: because current machines have memory limitations,
%            this only uses a randomly chosen subset of the data
%      Intended to be called by coaddm.m
% Outputs:
%     Only graphics, in current Figure

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



%  Set local parameters
%
ndecim = 350 ;   %  Number of randomly chosen points to reduced to
                 %  Set this to 0 to use all points (for better machines)
cdist = .3 ;  %  Parameter controlling random radii in randomly
              %  chosen grid of points, according to "paralellogram
              %  distribution":
              %  0 - "triangular distribution", i.e. weights to
              %          give unif dist'n on disk
              %  1 - "unif distribution", i.e. uniform on radius,
              %          so "many more near center
              %  between 0 and 1 -  a compromise
ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array
icut = 1 ;    %  1 to cut off first two elements in file read
ishowg = 0 ;     %  1 to show points used in interpolation
decseed = 29385534 ;    %  seed for random decimation



%  Set up global variables
%
%  First some general globals
global R0 NMAX ;
%  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

%  Next some globals for file names, etc.
global INFILESTR RFILESUF CFILESUF RISTR CISTR;
%  INFILESTR - Prefix of filename for input data
%  RFILESUF  - radii data file suffix 
%  CFILESUF  - curvature data file suffix
%  RISTR     - string to read in radius file
%  CISTR     - string to read in curvature file

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

%  Next the globals for handles to figure windows and plots
global DDMHAND RDHANDPOS ;
%  DDMHAND - handle for Figure Window with Display, Decomposition Menu
%  RDHANDPOS - Position for the new reconstructed plot


%  Read in data
%
% RISTR, CISTR (read-in file strings, set in advance in cornean.m)
fidr = fopen(RISTR) ;
fidc = fopen(CISTR) ;
          %  assigns a file handle, (read only assumed)
if (fidr == -1) | (fidc == -1) ;      %  then a file not found, so report

  disp('   ') ;
  disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
  disp('!!!   Can''t do raw data display, because can''t find file(s):') ;
  
  if (fidr == -1) ;
    disp(['!!!       ' RISTR]) ;
  end ;
 
  if (fidc == -1) ;
    disp(['!!!       ' CISTR]) ;
  end ;

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

else ;


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



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

  rdata = reshape(rdata,256,27) ;
  cdata = reshape(cdata,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) ;
  cdata = cdata(:,2:27) ;
          %  cutoff first columns of all zeros:





  %  Do random selection of elements from grid (if needed)
  %
  if ndecim ~= 0 ;
    rand('seed',decseed) ;
    randang = ceil(256 * rand(ndecim,1)) ;     
          %  Random index among 1,2,...,256
    %  hist(randang) ; pause ;
    %      Used this to test


      if cdist < 1 ;
        rpara = sqrt(cdist^2 + 4 * rand(ndecim,1) * (1 - cdist)) ;
        rpara = (-cdist + rpara) / (2 * (1 - cdist)) ;
      else ;
        rpara = rand(ndecim,1) ;
      end ;
    randrad = ceil(26 * rpara) ;
          %  Index among 1,2,...,26, with "triangular distribution",
          %  since sqrt(Unif) has c.d.f = x^2, i.e. density 2x
    %cdist
    %hist(randrad,1:26) ; pause ; 
    %      Used this to test

    vind = find(sparse(randang,randrad,ones(ndecim,1),256,26)) ;
          %  indices where have randomly chosen entries

      %  Turn into columns, decimate 
    rdata = rdata(:) ;
    rdata = rdata(vind) ;
    adata = adata(:) ;
    adata = adata(vind) ;
    cdata = cdata(:) ;
    cdata = cdata(vind) ;
    %disp(['After random decimation, have ' num2str(length(rdata)) ' values']) ;

  end ;




  %  remove inappropriate values
  %
  %  first missings in radii files
  vind = find(rdata == 0) ;
           %  vector of indices where have a missing
  rdata(vind) = [] ;
  adata(vind) = [] ;
  cdata(vind) = [] ;
           %  Deletes that part of the vectors

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

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




  %  Get x,y coordinates and grid the data
  %
  xdata = rdata .* cos(adata) ;
  ydata = rdata .* sin(adata) ;
  %disp('Now Grid the data') ;
    xgrid = linspace(-R0,R0,ngrid) ;
    ygrid = xgrid ;
            %  Vectors of image coordinate values 
  [mxgrid,mygrid,mout] = griddata(xdata,ydata,cdata,xgrid,ygrid') ;




  %  Setup Keratron curvature color map
  %
  %  These use the "full range" Keratron data, from the email message:
  %             cocolora.msg
  %   which were copied into the ASCII file 
  %             cocolor1.dat
  %  for importing into matlab programs
  %    locurve & hicurve are upper and lower ends of the "range"
  %    for each color
  eval(['load ' INDATPATH 'cocolor1.dat']) ;
          %  Assumes this file is in the directory with input data
  locurve = cocolor1(:,1) ;
  hicurve = cocolor1(:,2) ;
  cyan = cocolor1(:,3) ;
  magenta = cocolor1(:,4) ;
  yellow = cocolor1(:,5) ;
  black = cocolor1(:,6) ;

  red   = 1 - cyan + black ;
  green = 1 - magenta + black;
  blue  = 1 - yellow + black ;

  comap = [red green blue] ;
  ncomap = size(comap,1) ;
          %  number of rows in the color map ;

  %  make sure values are positive
  flag = comap < 0 ;
  fflag = find(flag) ;
  comap(fflag) = zeros(length(fflag),1) ;

  %  make sure values are less than 1
  flag = comap > 1 ;
  fflag = find(flag) ;
  comap(fflag) = ones(length(fflag),1) ;




  %  Create the matrix of image pixels
  %
  u01 = linspace(0,1,ncomap+1) ;
          %  integers representing the Keratron colors
  colcurv = interp1([-100; hicurve], u01, mout(:)) ;
          %  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)) ;  



  %  Crop to image only inside circle
  %
  mcrop = ones(ngrid) ;            
    vind = find(mxgrid.^2 + mygrid.^2 > R0^2) ;
                %  Vector of indices where image is outside circle

  mcrop(vind) = nan * ones(length(vind),1) ;
                %  Cropping Matrix, multiply by this to cut off images
                %  outside circle of radius r_0.       
  moutc = mout .* mcrop ;






  %  Create New Window for This plot
  %
  rdhand = figure ;
  set(rdhand,'Position', RDHANDPOS) ;




  %  Create Image
  %
  pcolor(xgrid,ygrid,moutc) ;
    shading flat ;
    axis('square') ;
      titstr = ['Raw Curvature in ' INFILESTR '.tla'] ;
    title(titstr) ;
      caxis([0 1]) ;
    if ishowg == 1 ;     %  Then add random grid have decimated to
      hold on ;
        polar(adata,rdata,'ok') ;
        axis('off') ;
      hold off ;
    end ;
    colormap(comap) ;


end ;



