disp('Running MATLAB script file lmtzcf12.m') ;
%
%    Makes Figure 1.2, for the Locantore, Marron, Tripoli,
%    Zhang, Cohen paper.
%


%  use lines from coarwad.m
%
%  First some general globals
global INDATPATH ;
%
%  First define variables (nothing is local)
R0 = 4.0 ;
NMAX = 10 ;
IBASIS = 1 ;
IRECENT = 0 ;
CRAD = 0 ;
CANG = 0 ;

INDATPATH = '\matlab\cornean\indata\' ;
SUBDIRSTR = 'milan\' ;
RFILESUF = '.xlb' ;
          %  radii data file suffix 
CFILESUF = '.tlb' ;
          %  curvature data file suffix


for ipart = 1:1 ;    %  loop through one image

  if ipart == 1 ;
    INFILESTR = '000121_2' ;   
  end ;


  RISTR = [INDATPATH SUBDIRSTR INFILESTR RFILESUF] ;
  CISTR = [INDATPATH SUBDIRSTR INFILESTR CFILESUF] ;
  %  Set up Keratron colormap
  %
  [COMAP, HICURVE] = coacmap ;


  %  Set local parameters
  %
  ndecim = 0 ;
  %ndecim = 200 ;   %  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
  nang  = 256 ;    %  number of semi-meridians read per radial ring
  okread = 0 ;	 %  1 to show successful file read
  ishowg = 0 ;     %  1 to show points used in interpolation
  decseed = 29385534 ;    %  seed for random decimation




  %  Read in data
  %
  % RISTR, CISTR (read-in file strings, set in advance in cornean.m)

  fidr = fopen(RISTR, 'r', 'l') ;
	%  IEEE format for binary files
  fidc = fopen(CISTR, 'r', 'l') ;
          %  assigns a file handle, (read only assumed)
  if (fidr == -1) | (fidc == -1) ;      %  then a file not found, so reformat
    if (upper(RFILESUF) == RFILESUF);
       RFILESUF = lower(RFILESUF);		% handles case sensitivity
       CFILESUF = lower(CFILESUF);		% .xlb vs. .XLB
    else				
      RFILESUF = upper(RFILESUF);
      CFILESUF = upper(CFILESUF);
    end;
    newristr = [INDATPATH SUBDIRSTR INFILESTR RFILESUF];
    newcistr = [INDATPATH SUBDIRSTR INFILESTR CFILESUF];
    fidr = fopen(newristr, 'r', 'l');
    fidc = fopen(newcistr, 'r', 'l');
	    %   try binary open with different suffix
  end;

  if (fidr == -1) | (fidc == -1) ;      %  then a file not found, so try ascii
    newrsuf = upper([RFILESUF(1:3) 'a']);
    newcsuf = upper([CFILESUF(1:3) 'a']);
    newristr = [INDATPATH SUBDIRSTR INFILESTR newrsuf];
    newcistr = [INDATPATH SUBDIRSTR INFILESTR newcsuf];
    fidr = fopen(newristr);
    fidc = fopen(newcistr);
	    %   new file handles for ascii format

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

    if (fidr == -1) | (fidc == -1) ;   % then file not found, so give error
      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 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
      cdata = fscanf(fidc,'%g') ;
      fclose(fidr) ;
      fclose(fidc) ;
     
      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) = [];
        cdata(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);

    cstatus = fseek(fidc, 2, -1);
    cdata = fread(fidc, 'float');
    fclose(fidc);

    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) = [];
      cdata(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) ;
    cdata = reshape(cdata,nang,nring) ;

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

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


    %  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) ;
    cdata(vor0ind) = zeros(nor0ind,1) ;

  %#############
  % !!! Note:  rinner and cinner draw an inside line for possible missing
  % 	     data on the interior (reinhard).
  % inedger and inedgec have been commented out, and may be used in the
  % future for other types of analyses.

    %  construct vector of data for "outer ring"
    %
    [rinner, rdvmvind] = coamdf(rdata, nring-1, nang) ;
          %  vector of indices before first missing values for rdata
    [cinner, cdvmvind] = coamdf(cdata, nring-1, nang) ;
          %  vector of indices before first missing values for hdata
    vmvind = min([rdvmvind'; cdvmvind'])' ;
          %  vector of indices before first missing values (exterior)
    vinner = max([rinner'; cinner'])' ;
          %  vector of indices after first missing values (interior)
    edgerdata = [] ;
    edgecdata = [] ;
  %  inedger = [];
  %  inedgec = [];
    for iang = 1:nang ;
  %    inedger   = [inedger; rdata(iang, vinner(iang, 1))] ;
  %    inedgec   = [inedgec; cdata(iang, vinner(iang, 1))] ;
      edgerdata = [edgerdata; rdata(iang,vmvind(iang,1))] ;
      edgecdata = [edgecdata; cdata(iang,vmvind(iang,1))] ;
    end ;
          %  take values where are inside R0, 
          %  and are before first non-missing
    edgeadata = adata(:,1) ;
          %  each angle

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

    %  Do random selection of elements from grid (if needed)
    %
    if ndecim ~= 0 ;
      rand('seed',decseed) ;
      randang = ceil(nang * rand(ndecim,1)) ;     
          %  Random index among 1,2,...,nang
      %  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((nring-1) * rpara) ;
          %  Index among 1,2,...,(nring-1), with "triangular distribution",
          %  since sqrt(Unif) has c.d.f = x^2, i.e. density 2x
      %cdist
      %hist(randrad,1:(nring-1)) ; pause ; 
      %      Used this to test



      vind = find(sparse(randang,randrad,ones(ndecim,1),nang,(nring-1))) ;
          %  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) ;



      %  Add on an outer ring, at R0, or at last non-missing value
      %  first delete overlaps
      videlete = [] ;
      for i = 1:nang ;
        flag = (rdata == edgerdata(i)) & (adata == edgeadata(i)) ;
          %  one where have an overlap
        if sum(flag) > 0 ;    %  then have some overlap,
                            %  so kick out of edge vectors
          videlete = [videlete; i] ;
        end ;
      end ;
      edgerdata(videlete) = [] ;
      edgeadata(videlete) = [] ;
      edgecdata(videlete) = [] ;
      %  next concatenate
      rdata = [rdata; edgerdata] ;
      adata = [adata; edgeadata] ;
      cdata = [cdata; edgecdata] ;
 


      %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']) ;




    %  Get x,y coordinates and grid the data
    %
  %  xdata = rdata .* cos(adata) ;
  % ydata = rdata .* sin(adata) ;
  [xdata, ydata] = pol2cart(adata, rdata);
    %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') ;


    ncomap = size(COMAP,1) ;
          %  number of rows in the color map ;



    %  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 Image
    %
    subplot(1,1,ipart) ;
    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 ;


end ;


  orient tall ;
  set(gcf,'PaperPosition',[3.25, 7.0, 2.0, 2.0]) ;
        %  3.25 inches on left and right, 2.0 interior, 8.5 total
        %  7.0 from bottom, 2.0 interior, 11 total
print -deps \matlab\steve\ps\lmtzcf12.eps ;


