%MATLAB Script File: coraw1.m
%
%   For the Cornea Study, this plots raw images
%
clear ;       %  To avoid memory problems
clf reset ;
format compact ;
running = 'MATLAB Script File:  coraw1.m'

autop = 1 ;   %  1 to create a color postscript file

ihires = 1 ;  %  Controls resolution of plot:
              %    0 - for lo res 64 x64  (fastest)
              %    1 - to interpolate to hi res 256 x 256  (slower)

iimag = 9 ;   %  0 - original test case
              %  1 - 84_1
              %  2 - 85_1
              %  3 - 122_2
              %  4 - 121_2
              %  5 - 55_1
              %  6 - 28_2
              %  7 - 130_1
              %  8 - 139_2
              %  9 - 6_1
              %  10 - 53_1
              %  11 - 43_1

isurf = 3 ;   %  1 - for height map
              %  2 - for height map - sphere
              %  3 - for curvature map
              %             (doesn't work for iimag = 0)

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

        

r0 = 4.2 ;    %  Usual radius of analysis
ndecim = 400 ;%  Number of points to randomly decimate to
ndecim = 350 ;%  Number of points to randomly decimate to

ngrid = 128 ; %  Construct surfaces as ngrid x ngrid array
ngrid = 64 ; %  Construct surfaces as ngrid x ngrid array
%ngrid = 20 ; %  Construct surfaces as ngrid x ngrid array

ninterp = 128 ;  %  Number of points to interpolate to for final plot

xgrid = linspace(-r0,r0,ngrid) ;
ygrid = xgrid ;
              %  Vectors of image coordinate values 

hrxg = linspace(-r0,r0,ninterp) ;
hryg = hrxg ;


    %  Create input strings
if iimag == 0 ;
  patstr = '1st Test' ;
  ristr = '000063_2.xpa' ;
  histr = '000063_2.ypa' ;
  cstep = 10 ;  %  Step size for color scale
  spherad = 8 ;  %  Radius of sphere to subtract
elseif iimag == 1 ;
  patstr = 'Image 84_1' ;
  ristr = '000084_1.xla' ;
  histr = '000084_1.zla' ;
  cistr = '000084_1.tla' ;
  pdist = .123 ;
  pang = 123 ;
  spherad = 7.941 ;  %  Radius of sphere to subtract
  cstep = 6 ;  %  Step size for color scale
%  spherad = 337.5 / 43.5 ;  %  Radius of sphere to subtract (Diotpers)
%  cstep = 7 ;  %  Step size for color scale
elseif iimag == 2 ;
  patstr = 'Image 85_1' ;
  ristr = '000085_1.xla' ;
  histr = '000085_1.zla' ;
  cistr = '000085_1.tla' ;
  pdist = .254 ;
  pang = 68 ;
  spherad = 7.258 ; %  Radius of sphere to subtract
  cstep = 25 ;  %  Step size for color scale
elseif iimag == 3 ;
  patstr = 'Image 122_2' ;
  ristr = '000122_2.xla' ;
  histr = '000122_2.zla' ;
  cistr = '000122_2.tla' ;
  pdist = .686 ;
  pang = 89 ;
  spherad = 7.8488 ; %  Radius of sphere to subtract
  cstep = 35 ;  %  Step size for color scale
elseif iimag == 4 ;
  patstr = 'Image 121_2' ;
  ristr = '000121_2.xla' ;
  histr = '000121_2.zla' ;
  cistr = '000121_2.tla' ;
  pdist = .525 ;
  pang = 108 ;
  spherad = 337.5 / 41 ;  %  Radius of sphere to subtract (Diotpers)
  cstep = 40 ;  %  Step size for color scale
elseif iimag == 5 ;
  patstr = 'Image 55_1' ;
  ristr = '000055_1.xla' ;
  histr = '000055_1.zla' ;
  cistr = '000055_1.tla' ;
  pdist = .439 ;
  pang = 121 ;
  spherad = 337.5 / 43.5 ;  %  Radius of sphere to subtract (Diotpers)
  cstep = 12 ;  %  Step size for color scale
elseif iimag == 6 ;
  patstr = 'Image 28_2' ;
  ristr = '000028_2.xla' ;
  histr = '000028_2.zla' ;
  cistr = '000028_2.tla' ;
  pdist = .426 ;
  pang = 115 ;
  spherad = 8.0357 ; %  Radius of sphere to subtract
  cstep = 5 ;  %  Step size for color scale
elseif iimag == 7 ;
  patstr = 'Image 130_1' ;
  ristr = '000130_1.xla' ;
  histr = '000130_1.zla' ;
  cistr = '000130_1.tla' ;
  pdist = .458 ;
  pang = 112 ;
  spherad = 337.5 / 39.7 ;  %  Radius of sphere to subtract (Diotpers)
  cstep = 10 ;  %  Step size for color scale
%  spherad = 337.5 / 37.9 ;  %  Radius of sphere to subtract (Diotpers)
elseif iimag == 8 ;
  patstr = 'Image 139_2' ;
  ristr = '000139_2.xla' ;
  histr = '000139_2.zla' ;
  cistr = '000139_2.tla' ;
  pdist = .417 ;
  pang = 114 ;
  spherad = 8.4375 ; %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 9 ;
  patstr = 'Image 6_1' ;
  ristr = '000006_1.xla' ;
  histr = '000006_1.zla' ;
  cistr = '000006_1.tla' ;
  pdist = .454 ;
  pang = 116 ;
  spherad = 7.7586 ; %  Radius of sphere to subtract
  cstep = 6 ;  %  Step size for color scale
elseif iimag == 10 ;
  patstr = 'Image 53_1' ;
  ristr = '000053_1.xla' ;
  histr = '000053_1.zla' ;
  cistr = '000053_1.tla' ;
  pdist = .375 ;
  pang = 113 ;
  spherad = 8.766 ; %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 11 ;
  patstr = 'Image 43_1' ;
  ristr = '000043_1.xla' ;
  histr = '000043_1.zla' ;
  cistr = '000043_1.tla' ;
  pdist = .557 ;
  pang = 91 ;
  spherad = 8.5443 ; %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
end ;


pangr = (2 * pi * pang) / 360 ;
          %  Convert to radians


    %  Read in radii of data.
fid = fopen(['\matlab\steve\cornea\rawdata\' ristr]) ;
              %  opens file for reading
rdata = fscanf(fid,'%g') ;
              %  formatted read of file, assumes single column
if iimag >= 1 ;
              %  then need to cut off first two elements
  rdata = rdata(3:length(rdata)) ;
end ;
rdata = reshape(rdata,256,27) ;
              %  put in matrix form
fclose(fid) ;

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

      %  Read in rest of data.
if isurf == 1 | isurf == 2 ;     %  Then work with height data
  surfstr = ' Heights' ;
  fid = fopen(['\matlab\steve\cornea\rawdata\' histr]) ;
                %  opens file for reading
  hdata = fscanf(fid,'%g') ;
                %  formatted read of file, assumes single column
  if iimag >= 1 ;
                %  then need to cut off first two elements
    hdata = hdata(3:length(hdata)) ;
  end ;
  hdata = reshape(hdata,256,27) ;
                %  put in matrix form
  fclose(fid) ;

  if isurf == 2 ;   %  Then need to subtract sphere
    surfstr = [surfstr ' - Sphere'] ;
    sphht = spherad - sqrt(spherad.^2 - rdata.^2) ;
    hdata = hdata + sphht ;
  end ;

elseif isurf == 3 ;     %  Then work with curvature data
  surfstr = ' Curvatures' ;

  fid = fopen(['\matlab\steve\cornea\rawdata\' cistr]) ;
                %  opens file for reading
  hdata = fscanf(fid,'%g') ;
                %  formatted read of file, assumes single column
  if iimag >= 1 ;
                %  then need to cut off first two elements
    hdata = hdata(3:length(hdata)) ;
  end ;
  hdata = reshape(hdata,256,27) ;
                %  put in matrix form
  fclose(fid) ;

end ;

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

%for iang = 



  %  Do random selection of elements from grid
rand('seed',29385534) ;
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) ;
hdata = hdata(:) ;
hdata = hdata(vind) ;
disp(['After random decimation, have ' num2str(length(rdata)) ' values']) ;

if isurf <= 2 ;
          %  Then add 0 point back in
  rdata = [0; rdata] ;
  adata = [0; adata] ;
  hdata = [0; hdata] ;
end ;

    %  remove the missing values in radii files
vind = find(rdata == 0) ;
             %  vector of indices where have a missing
if isurf <= 2 ;
             %  if 0 added back in
  vind = vind(2:length(vind)) ;
end ;
rdata(vind) = [] ;
adata(vind) = [] ;
hdata(vind) = [] ;
             %  Deletes that part of the vectors

    %  remove the missing values in curvature files
vind = find(hdata == 0) ;
             %  vector of indices where have a missing
if isurf <= 2 ;
             %  if 0 added back in
  vind = vind(2:length(vind)) ;
end ;
rdata(vind) = [] ;
adata(vind) = [] ;
hdata(vind) = [] ;
             %  Deletes that part of the vectors

disp(['After removing missings, have ' num2str(length(rdata)) ' values']) ;


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


    %  get x,y coordinates
xdata = rdata .* cos(adata) ;
ydata = rdata .* sin(adata) ;


disp('Now Grid the data') ;
[mxgrid,mygrid,mout] = griddata(xdata,ydata,hdata,xgrid,ygrid') ;



if isurf == 2 ;
          %  Then setup Keratron height 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
  load cornea\cocolor1.dat ;
  cyan = cocolor1(:,3) ;
  magenta = cocolor1(:,4) ;
  yellow = cocolor1(:,5) ;
  black = cocolor1(:,6) ;
  %  Don't use columns 1 and 2, because those are for curvature

  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
  %     (since twiddling "black" can leave this range)
  flag = comap < 0 ;
  fflag = find(flag) ;
  comap(fflag) = zeros(length(fflag),1) ;

  %  make sure values are less than 1
  %     (since twiddling "black" can leave this range)
  flag = comap > 1 ;
  fflag = find(flag) ;
  comap(fflag) = ones(length(fflag),1) ;

  colormap(comap) ;

      cstep = cstep / 1000 ;
          %  cstep is "step size in color scale", different for
          %  each image

      u01 = linspace(0,1,ncomap+1) ;
        lbd = -11.5 * cstep ;
        ubd = lbd + ncomap * cstep ;
      colgrid = linspace(lbd,ubd,ncomap+1) ;
          %  equally spaced grid on desired color scale

    vmout = mout(:) ;
          %  vector version for easier manipulation
          %  NOTE: interp1 seems to need vector argument
      indlo = find(vmout < lbd) ;
    vmout(indlo) = lbd * ones(length(indlo),1) ;
      indup = find(vmout > ubd) ;
    vmout(indup) = lbd * ones(length(indup),1) ;
          %  NOTE: interp1 seems to need values within grid

    colcurv = interp1(colgrid, 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.
  mout = reshape(colcurv,size(mout)) ;  

elseif isurf == 3 ;
          %  Then 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
  load cornea\cocolor1.dat ;
  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) ;

  colormap(comap) ;

    u01 = linspace(0,1,ncomap+1) ;
    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)) ;  
end ;



disp('Now crop and make plot') ;
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 ;



if ihires == 0 ;
      %  Make low resolution plot
  pcolor(mxgrid,mygrid,moutc) ;
    shading flat ;
    axis('square') ;
    title([patstr surfstr]) ;
    if isurf == 2 ;   %  if heightmap - sphere
      caxis([0 1]) ;
    elseif isurf == 3 ;   %  if curvature map
  %    caxis([14 55]) ;    %  old color scheme
      caxis([0 1]) ;
    end ;
    hold on ;
      polar(adata,rdata,'ok') ;
            %  Add random grid have decimated too
      polar((pi + pangr),pdist,'*k') ;
            %  Show Keratron center
      axis('off') ;
    hold off ;

elseif ihires == 1 ;
      %  Make high resolution plot
  [mhrxg,mhryg] = meshgrid(hrxg,hryg) ;
  mhrout = interp2(mxgrid,mygrid,mout,mhrxg,mhryg) ;
  mcrop = ones(ninterp) ;            
    vind = find(mhrxg.^2 + mhryg.^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.       
  mhrout = mhrout .* mcrop ;
  
      %  Make high resolution plot
  pcolor(mhrxg,mhryg,mhrout) ;
    shading flat ;
    axis('square') ;
    title([patstr surfstr]) ;
    if isurf == 3 ;   %  if curvature map
  %    caxis([14 55]) ;    %  old color scheme
      caxis([0 1]) ;
    end ;
    hold on ;
      polar(adata,rdata,'ok') ;
            %  Add random grid have decimated too
      polar((pi + pangr),pdist,'*k') ;
            %  Show Keratron center
      axis('off') ;
    hold off ;
end ;



if autop == 1 ;
  if isurf == 2 ;
    eval(['print \matlab\steve\ps\coraw' num2str(iimag) 'h.ps -dpsc']) ;
  elseif isurf == 3 ;
    eval(['print \matlab\steve\ps\coraw' num2str(iimag) 'c.ps -dpsc']) ;
  end ;
end ;


