disp('MATLAB Script File:  cosd1.m') ;
%
%   For the Cornea Study, this does a "surfaces as data" approach
%       to understanding the variability in the corneal groups
%
format compact ;

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

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

isurf = 3 ;   %  1 - for height map
              %  2 - for height map - sphere
              %  3 - for curvature map

ibas = 1 ;    %  1 - Radial Fourier - Legendre
              %  2 - Diameter Fourier - Legendre
              %            (not working yet)
              %  3 - Zernike
              %            (not working yet)


r0 = 4.2 ;    %  Usual radius of analysis

ngrid = 64 ; %  Construct surfaces as ngrid x ngrid array
if ihires == 2 ;
          %  then change ngrid to larger value
  ngrid = 256 ; %  Construct surfaces as ngrid x ngrid array
end ;
hrngrid = 256 ; %  Interpolate to hrngrid x hrngrid for 
                %  higher resolution part

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) ;
              %  Cropping Matrix, multiply by this to cut off images
              %  outside circle of radius r_0.       

nr = 10 ;
nc = 7 ;

if igroup == 1 ;
  namestr = '\matlab\steve\cornea\fm\fname1mo.asc' ;
  imagstr = '1 Mo.' ;
elseif igroup == 2 ;
  namestr = '\matlab\steve\cornea\fm\fnameker.asc' ;
  imagstr = 'Ker.' ;
end ;


disp('Loading feature matrices for this group') ;
fid = fopen(namestr,'r') ;
vfindx = fscanf(fid,'%4d') ;

prestr = '\matlab\steve\cornea\fm\fm42' ;
mfv = [] ;
for i = 1:size(vfindx,1) ;
  if vfindx(i) < 10 ;
    eval(['load ' prestr '000' num2str(vfindx(i))]) ;
  elseif vfindx(i) < 100 ;
    eval(['load ' prestr '00' num2str(vfindx(i))]) ;
  elseif vfindx(i) < 1000 ;
    eval(['load ' prestr '0' num2str(vfindx(i))]) ;
  else ;
    eval(['load ' prestr num2str(vfindx(i))]) ;
  end ;
          %  This creates a number of variables, especially fm,
          %  a 19 x 7 "feature matrix"
  mfv = [mfv fm(:)] ;
          %  put matrix into a column vector and add to right
end ;

disp(['feature matrix has size ' num2str(size(mfv,1)) ' rows ']) ;
disp(['and ' num2str(size(mfv,2)) ' cols']) ;


%%  These lines tested that these are the correct feature matrices.
%%  First load data from the file a42p0592.tst, which is an edited version
%%  of a42p0592.dat, to remove headers, and put all matrix rows on the 
%%  same line.
%load a42p0592.tst ;
%testamp = a42p0592(1:10,:) ;
%testpha = a42p0592(11:20,:) * pi / 180 ;
%testfm = testamp(1,:) .* cos(testpha(1,:)) ;
%for i = 2:10 ;
%  testfm = [testfm; testamp(i,:) .* cos(testpha(i,:))] ;
%  testfm = [testfm; testamp(i,:) .* sin(testpha(i,:))] ;
%end ;
%testfm = testfm(:) ;
%disp('Check this is essentially 0') ;
%max(abs(testfm - mfv(:,1))) 
%pause


disp('Working on eigenector decomposition') ;
d = size(mfv,1) ;
n = size(mfv,2) ;
sst = sum(sum(mfv.^2)) ;
%  get mean feature vector, and resids from that
vmean = mean(mfv') ;
vsd = std(mfv') ;
ssm = n * sum(vmean.^2) ;
ssmr = sst - ssm ;
mresid = mfv - vec2mat(vmean',n) ;

%  Get cov's and eigen values
mcov = cov(mresid') ;
          %  Get covariance matrix, transpose, since want "coordinates
          %                      as variables"

[meigvec,meigval] = eig(mcov) ;
veigval = (diag(meigval)) ;
          %  Get eigenvalues and eigenvectors

[veigval,vind] = sort(veigval) ;
meigvec = meigvec(:,vind) ;
          %  Sort to avoid possible ordering problems
  
ev1 = meigvec(:,d) ;
          %  eigenvector with largest eigenvalue
pc1 = mresid' * ev1 ;
          %  vector of projections of data onto ev1
pc1rep = median(abs(pc1 - mean(pc1))) ;
          %  "most representative pc", defined as MAD
mpc1 = ev1 * pc1rep ;
          %  "average vector" for display
mresid1 = (mresid - ev1 * pc1') ;
ss1r = sum(sum(mresid1.^2)) ;
ss1 = ssmr - ss1r ;

ev2 = meigvec(:,d-1) ;
          %  eigenvector with largest eigenvalue
pc2 = mresid1' * ev2 ;
          %  vector of projections of data onto ev2
pc2rep = median(abs(pc2 - mean(pc2))) ;
          %  "most representative pc", defined as MAD
mpc2 = ev2 * pc2rep ;
          %  "average vector" for display
mresid2 = (mresid1 - ev2 * pc2') ;
ss2r = sum(sum(mresid2.^2)) ;
ss2 = ss1r - ss2r ;

ev3 = meigvec(:,d-2) ;
          %  eigenvector with largest eigenvalue
pc3 = mresid2' * ev3 ;
          %  vector of projections of data onto ev3
pc3rep = median(abs(pc3 - mean(pc3))) ;
          %  "most representative pc", defined as MAD
mpc3 = ev3 * pc3rep ;
          %  "average vector" for display
mresid3 = (mresid2 - ev3 * pc3') ;
ss3r = sum(sum(mresid3.^2)) ;
ss3 = ss2r - ss3r ;

ev4 = meigvec(:,d-3) ;
          %  eigenvector with largest eigenvalue
pc4 = mresid3' * ev4 ;
          %  vector of projections of data onto ev4
pc4rep = median(abs(pc4 - mean(pc4))) ;
          %  "most representative pc", defined as MAD
mpc4 = ev4 * pc4rep ;
          %  "average vector" for display
mresid4 = (mresid3 - ev4 * pc4') ;
ss4r = sum(sum(mresid4.^2)) ;
ss4 = ss3r - ss4r ;


fclose('all') ;


%  Reshape representative guys to feature matrix form
fmmean = reshape(vmean,19,7) ;
fmpc1 = reshape(mpc1,19,7) ;
fmpc2 = reshape(mpc2,19,7) ;
fmpc3 = reshape(mpc3,19,7) ;
fmpc4 = reshape(mpc4,19,7) ;


%  start plotting stuff, from corecon5.m
disp('Generating surface') ; 
if isurf == 1 | isurf == 2 ;   %  then work with height map
  surfstr = ' Heights' ;

  if ibas == 1 ;     %  then use the radial Fourier - Legendre basis

    moutm = corh(fmmean,mrgrid,mthetag,r0) ;
    mout1 = corh(fmpc1,mrgrid,mthetag,r0) ;
    mout2 = corh(fmpc2,mrgrid,mthetag,r0) ;
    mout3 = corh(fmpc3,mrgrid,mthetag,r0) ;
    mout4 = corh(fmpc4,mrgrid,mthetag,r0) ;
          %  Radial Fourier - Legendre reconstructions of Heights

  elseif ibas == 2 ;     %  then use the diameter Fourier - Legendre basis

  elseif ibas == 3 ;     %  then use the Zernike basis

  end ;

  if isurf == 2 ;   %  Then need to subtract sphere
    surfstr = [surfstr ' - Sphere'] ;

%    minmout = min(min(mout))
%    maxmout = max(max(mout))
            %  use these lines to decide which sphere to subtract

    sphht = spherad - sqrt(spherad.^2 - mrgrid.^2) ;

%    minsph = min(min(sphht))
%    maxsph = max(max(sphht))
%    pause
            %  use these lines to decide which sphere to subtract


    moutm = moutm + sphht ;
    mout1 = mout1 + sphht ;
    mout2 = mout2 + sphht ;
    mout3 = mout3 + sphht ;
    mout4 = mout4 + sphht ;
            %  careful about whether this is + or -
  end ;

elseif isurf == 3 ;   %  then work with curvature map
  surfstr = ' Curvatures' ;

  if ibas == 1 ;     %  then use the radial Fourier - Legendre basis

    [mhm,mhpm,mhppm] = corhpp(fmmean',mrgrid,mthetag,r0) ;
    [mh1,mhp1,mhpp1] = corhpp(fmpc1',mrgrid,mthetag,r0) ;
    [mh2,mhp2,mhpp2] = corhpp(fmpc2',mrgrid,mthetag,r0) ;
    [mh3,mhp3,mhpp3] = corhpp(fmpc3',mrgrid,mthetag,r0) ;
    [mh4,mhp4,mhpp4] = corhpp(fmpc4',mrgrid,mthetag,r0) ;

  elseif ibas == 2 ;     %  then use the diameter Fourier - Legendre basis

  elseif ibas == 3 ;     %  then use the Zernike basis
  end ;

  moutm = -mhppm ./ (1 + mhpm.^2).^(3/2) ;
  mout1 = -mhpp1 ./ (1 + mhp1.^2).^(3/2) ;
  mout2 = -mhpp2 ./ (1 + mhp2.^2).^(3/2) ;
  mout3 = -mhpp3 ./ (1 + mhp3.^2).^(3/2) ;
  mout4 = -mhpp4 ./ (1 + mhp4.^2).^(3/2) ;
  moutm = moutm * 337.5 ;
  mout1 = mout1 * 337.5 ;
  mout2 = mout2 * 337.5 ;
  mout3 = mout3 * 337.5 ;
  mout4 = mout4 * 337.5 ;
          %  put on scale of diopters

end ;




%  do explicit flip of x's and y's, since pcolor wants:
%            rows as y's,   cols as x's
moutm = moutm' ;
mout1 = mout1' ;
mout2 = mout2' ;
mout3 = mout3' ;
mout4 = mout4' ;




if ihires == 1 ;
          %  Then compute higher resolution stuff
  hrxg = linspace(-r0,r0,hrngrid) ;
  hryg = hrxg ;
  [mhrxg,mhryg] = meshgrid(hrxg,hryg) ;
  mhroutm = interp2(mxgrid,mygrid,moutm,mhrxg,mhryg) ;
  mhrout1 = interp2(mxgrid,mygrid,mout1,mhrxg,mhryg) ;
  mhrout2 = interp2(mxgrid,mygrid,mout2,mhrxg,mhryg) ;
  mhrout3 = interp2(mxgrid,mygrid,mout3,mhrxg,mhryg) ;
  mhrout4 = interp2(mxgrid,mygrid,mout4,mhrxg,mhryg) ;
  mcrop = ones(hrngrid) ;            
  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.       
  mhroutm = mhroutm .* mcrop ;
  mhrout1 = mhrout1 .* mcrop ;
  mhrout2 = mhrout2 .* mcrop ;
  mhrout3 = mhrout3 .* mcrop ;
  mhrout4 = mhrout4 .* mcrop ;
end ;






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.
  moutm = reshape(colcurv,size(moutm)) ;  
  mout1 = reshape(colcurv,size(mout1)) ;  
  mout2 = reshape(colcurv,size(mout2)) ;  
  mout3 = reshape(colcurv,size(mout3)) ;  
  mout4 = reshape(colcurv,size(mout4)) ;  

elseif isurf == 3 ;
          %  Then setup Keratron 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) ;

    vmoutm = moutm(:) ;
    vmout1 = mout1(:) ;
    vmout1 = 45 * vmout1 / max(abs(vmout1)) + 43.75 ;
            %  makes largest magnitude = 45, then shifts 0 to 
            %  43.75, which is yellow color
    vmout2 = mout2(:) ;
    vmout2 = 45 * vmout2 / max(abs(vmout2)) + 43.75 ;
    vmout3 = mout3(:) ;
    vmout3 = 45 * vmout3 / max(abs(vmout3)) + 43.75 ;
    vmout4 = mout4(:) ;
    vmout4 = 45 * vmout4 / max(abs(vmout4)) + 43.75 ;
          %  vector version for easier manipulation
          %  NOTE: interp1 seems to need vector argument
      indlo = find(vmoutm < 0) ;
    vmoutm(indlo) = 0 * ones(length(indlo),1) ;
      indup = find(vmoutm > max(hicurve)) ;
    vmoutm(indup) = max(hicurve) * ones(length(indup),1) ;
          %  NOTE: interp1 seems to need values within grid
      indlo = find(vmout1 < 0) ;
    vmout1(indlo) = 0 * ones(length(indlo),1) ;
      indup = find(vmout1 > max(hicurve)) ;
    vmout1(indup) = max(hicurve) * ones(length(indup),1) ;
      indlo = find(vmout2 < 0) ;
    vmout2(indlo) = 0 * ones(length(indlo),1) ;
      indup = find(vmout2 > max(hicurve)) ;
    vmout2(indup) = max(hicurve) * ones(length(indup),1) ;
      indlo = find(vmout3 < 0) ;
    vmout3(indlo) = 0 * ones(length(indlo),1) ;
      indup = find(vmout3 > max(hicurve)) ;
    vmout3(indup) = max(hicurve) * ones(length(indup),1) ;
      indlo = find(vmout4 < 0) ;
    vmout4(indlo) = 0 * ones(length(indlo),1) ;
      indup = find(vmout4 > max(hicurve)) ;
    vmout4(indup) = max(hicurve) * ones(length(indup),1) ;

  colormap(comap) ;

    u01 = linspace(0,1,ncomap+1) ;
    colcurv = interp1([0; hicurve], u01, vmoutm) ;
          %  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
  moutm = reshape(colcurv,size(moutm)) ;  
    colcurv = interp1([0; hicurve], u01, vmout1) ;
  mout1 = reshape(colcurv,size(mout1)) ;  
    colcurv = interp1([0; hicurve], u01, vmout2) ;
  mout2 = reshape(colcurv,size(mout2)) ;  
    colcurv = interp1([0; hicurve], u01, vmout3) ;
  mout3 = reshape(colcurv,size(mout3)) ;  
    colcurv = interp1([0; hicurve], u01, vmout4) ;
  mout4 = reshape(colcurv,size(mout4)) ;  
end ;




if ihires == 1 ;
          %  Then compute higher resolution stuff
  hrxg = linspace(-r0,r0,hrngrid) ;
  hryg = hrxg ;
  [mhrxg,mhryg] = meshgrid(hrxg,hryg) ;
  mhroutm = interp2(mxgrid,mygrid,moutm,mhrxg,mhryg) ;
  mhrout1 = interp2(mxgrid,mygrid,mout1,mhrxg,mhryg) ;
  mhrout2 = interp2(mxgrid,mygrid,mout2,mhrxg,mhryg) ;
  mhrout3 = interp2(mxgrid,mygrid,mout3,mhrxg,mhryg) ;
  mhrout4 = interp2(mxgrid,mygrid,mout4,mhrxg,mhryg) ;
  mcrop = ones(hrngrid) ;            
  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.       
  mhroutm = mhroutm .* mcrop ;
  mhrout1 = mhrout1 .* mcrop ;
  mhrout2 = mhrout2 .* mcrop ;
  mhrout3 = mhrout3 .* mcrop ;
  mhrout4 = mhrout4 .* mcrop ;
else ;
  moutm = moutm .* mcrop ;
          %  Cut off part outside analysis circle
  mout1 = mout1 .* mcrop ;
  mout2 = mout2 .* mcrop ;
  mout3 = mout3 .* mcrop ;
  mout4 = mout4 .* mcrop ;
end ;




disp('Making Mean Plot') ;
subplot(3,3,1) ;

if ihires == 1 ;
  pcolor(mhrxg,mhryg,mhroutm) ;
else ;
  pcolor(mxgrid,mygrid,moutm) ;
end ;
  shading flat ;
  axis('square') ;
  title([imagstr surfstr ' mean']) ;
    ssstr = ['% of SS,T = ' num2str(100 * ssm / sst)] ;
  text(-2.4,3.6,ssstr) ;  
  if isurf == 2 | isurf == 3 ;   
          %  if heightmap - sphere or curvature
    caxis([0 1]) ;
  end ;




disp('Working on 1st eigenvector') ;
subplot(3,3,2) ;

if ihires == 1 ;
  pcolor(mhrxg,mhryg,mhrout1) ;
else ;
  pcolor(mxgrid,mygrid,mout1) ;
end ;
  shading flat ;
  axis('square') ;
  title([imagstr surfstr ' pc1']) ;
    ssstr = ['% of SS,MR = ' num2str(100 * ss1 / ssmr)] ;
  text(-2.4,3.6,ssstr) ;  
%  if isurf == 2 | isurf == 3 ;   
%          %  if heightmap - sphere or curvature
%    caxis([0 1]) ;
%  end ;

%  plot PC1 Components
subplot(3,3,3) ;
[kde, xgrid] = gpkde(pc1) ;
plot(xgrid,kde,'w') ;
  vlim = axis ;
          %  Get axis values for drawing vertical lines
  bottom = vlim(3) ;
  top = vlim(4) ;
hold on ;
  plot([pc1 pc1], [bottom top],'b') ;
  plot(xgrid,kde,'w') ;
          %  Do a replot to ensure seeing kde
  title('Projections, PC1') ;
hold off ;





disp('Working on 2nd eigenvector') ;
subplot(3,3,5) ;

if ihires == 1 ;
  pcolor(mhrxg,mhryg,mhrout2) ;
else ;
  pcolor(mxgrid,mygrid,mout2) ;
end ;
  shading flat ;
  axis('square') ;
  title([imagstr surfstr ' pc2']) ;
    ssstr = ['% of SS,1R = ' num2str(100 * ss2 / ss1r)] ;
  text(-2.4,3.6,ssstr) ;  
%  if isurf == 2 | isurf == 3 ;   
%          %  if heightmap - sphere or curvature
%    caxis([0 1]) ;
%  end ;

%  plot PC2 Components
subplot(3,3,6) ;
[kde, xgrid] = gpkde(pc2) ;
plot(xgrid,kde,'w') ;
  vlim = axis ;
          %  Get axis values for drawing vertical lines
  bottom = vlim(3) ;
  top = vlim(4) ;
hold on ;
  plot([pc2 pc2], [bottom top],'b') ;
  plot(xgrid,kde,'w') ;
          %  Do a replot to ensure seeing kde
  title('Projections, PC2') ;
hold off ;




disp('Working on 3rd eigenvector') ;
subplot(3,3,8) ;

if ihires == 1 ;
  pcolor(mhrxg,mhryg,mhrout3) ;
else ;
  pcolor(mxgrid,mygrid,mout3) ;
end ;
  shading flat ;
  axis('square') ;
  title([imagstr surfstr ' pc3']) ;
    ssstr = ['% of SS,2R = ' num2str(100 * ss3 / ss2r)] ;
  text(-2.4,3.6,ssstr) ;  
%  if isurf == 2 | isurf == 3 ;   
%          %  if heightmap - sphere or curvature
%    caxis([0 1]) ;
%  end ;

%  plot PC3 Components
subplot(3,3,9) ;
[kde, xgrid] = gpkde(pc3) ;
plot(xgrid,kde,'w') ;
  vlim = axis ;
          %  Get axis values for drawing vertical lines
  bottom = vlim(3) ;
  top = vlim(4) ;
hold on ;
  plot([pc3 pc3], [bottom top],'b') ;
  plot(xgrid,kde,'w') ;
          %  Do a replot to ensure seeing kde
  title('Projections, PC3') ;
hold off ;




disp('Working on 4th eigenvector') ;
subplot(3,3,4) ;

if ihires == 1 ;
  pcolor(mhrxg,mhryg,mhrout4) ;
else ;
  pcolor(mxgrid,mygrid,mout4) ;
end ;
  shading flat ;
  axis('square') ;
  title([imagstr surfstr ' pc4']) ;
    ssstr = ['% of SS,3R = ' num2str(100 * ss4 / ss3r)] ;
  text(-2.4,3.6,ssstr) ;  
%  if isurf == 2 | isurf == 3 ;   
%          %  if heightmap - sphere or curvature
%    caxis([0 1]) ;
%  end ;

%  plot PC4 Components
subplot(3,3,7) ;
[kde, xgrid] = gpkde(pc4) ;
plot(xgrid,kde,'w') ;
  vlim = axis ;
          %  Get axis values for drawing vertical lines
  bottom = vlim(3) ;
  top = vlim(4) ;
hold on ;
  plot([pc4 pc4], [bottom top],'b') ;
  plot(xgrid,kde,'w') ;
          %  Do a replot to ensure seeing kde
  title('Projections, PC4') ;
hold off ;



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



