function sss1(data,vitype,varp,alpha,vhgp,vgp,bdryp) 
% SSS1, Significance in Scale Space, for finding features in 2d
%     Uses the best methods found in SiZer/2d study
%     based on slopes, curvatures or both
%
%     For each bandwidth, creates a gray level map of smoothed data
%     and overlays matrix of symbols, showing significance
%
%     This is based on circular (i.e. same bandwidth in all directions)
%     Gaussian kernel convolution smoothing
%
%   Can use first 1, 2, 3, 4, 5, 6, 7 or 8 arguments.
% Inputs:
%     data   - either n x m matrix of image data  
%                            (requires m > 2, otherwise interpreted as 
%                             density estimation data)
%                  or n x 2 matrix of 2d density estimation data
%                            X's in first column,  Y's in second
%    vitype  - vector controlling type of plot
%                  1st entry:
%                       0 - just do smooths only (same format as others)
%                       1 - arrows based on gradient only (default, when this
%                                                           is not specified)
%                       2 - dots for curvature only
%                       3 - combine arrows and dots (show dot only when no
%                                                     significant arrow)
%                       4 - Streamlines, showing only gradient information
%                  2nd entry:
%                       1 - use 1x1 grid (i.e. single pixels)
%                       2 - use 2x2 boxes (bigger symbols, default, when this
%                                                           is not specified)
%                                 (this has no effect when only smoothing or
%                                                          using streamlines)
%                  3nd entry:
%                       1 - static plot, 2x2 array (for 4 bandwidths)
%                                     (default, when this is not specified)
%                                      or single image (esp. for large images)
%                       2 - movie, looping through bandwidths
%                       3 - plot with slider (allows choice to print)
%                  4th entry:
%                       0 - don't display messages showing progress
%                       1 - display progress messages (default when this 
%                                                           is not specified)
%                  5th entry:
%                       0 - give warning and option prompt, when image size 
%                                is too large (default)
%                       1 - just continue with 2 x 2 array of plots,
%                                regardless of image size, no prompts
%                       2 - automatically switch to single page plots,
%                                when image size is too large for 2 x 2 array.
%                  6th entry:
%                       0 - draw circles for small ESS in regression,
%                                      but not for density estimation
%                                      (default)
%                       1 - no circles for small ESS in regression,
%                                      but draw for density estimation
%    varp    -  variance handling parameter
%                  for density estimation, this has no effect 
%                       [0,sig2] - use input known value of noise, sigma^2,
%                       [1,0]    - use pooled estimate of noise  (default)
%                       [1,1]    - use local estimate of noise
%    alpha   -  Usual "level of significance".  I.e. C.I.s have coverage
%                  probability 1 - alpha.  (0.05 when not specified)
%     vhgp   - vector of h (bandwidth) grid parameters:
%                       0 (or not specified)  -  use range:
%                                                 1 - 8 for images
%                                                 2 - 16 for densities
%                                   and 4 h's (for static plot)
%                                       25 (for movies and slider)
%                       > 0  (scalar)  -  Use just that h
%                       [hmin; hmax]  -  use hmin to hmax
%                                   and 4 h's (for static plot)
%                                       25 (for movies and slider)
%                 (Note: h less than 1 gives too small Eff. Sample Size)
%                       [hmin; hmax; nh]  -  use hmin to hmax and nh 
%                                    bandwidths (but no effect for static plot)
%     vgp    - vector of grid parameters:
%                  for image data, this has no effect (since size of
%                            data determines this, and use pixel numbers)
%                  for density estimation:
%                       0 (or not specified)  -  use maxes and mins
%                                of data and 64 x 64 grid, and linear binning
%                       1  -  use maxes and mins of data and 64 x 64 grid, 
%                                and nearest neighbor binning
%                       [xmin, xmax, nxg, ymin, ymax, nyg]  - use these 
%                                input values and linear binning
%                       [xmin, xmax, nxg, ymin, ymax, nyg, 1]  - use these 
%                                input values and nearest neighbor binning
%    bdryp   - boundary parameter:
%                  for image data:
%                       0  -  No boundary adjustment
%                       1 (or not specified) -  Adjust for boundary effects,
%                                     by subtracting the overall mean
%                       2  -  Adjust for boundary effects,
%                                     by subtracting the overall median
%                  for density estimation:
%                       0  -  truncate data outside range
%                       1 (or not specified)  -  move data outside range to
%                                   nearest endpoint
%
%  Not input as an argument, but CAN be used a global variable:
%
%  GRAPHHAND - Global Graphics Handle:
%                  If this is a global figure handle, then graphics will
%                       go to that figure, 
%                  This can be a vector of figure handles, so slider plots
%                       will not destroy previous sliders (if the next
%                       handle is appended)
%                  If this is an axes handle, and nh = 1, static graphics
%                       will be put in that axes, not the current figure
%          Caution:  If using this gives gives unpredictable results, try:
%                      close all ;    clear all ;
%
% Output:
%     Draws a 2x2 array of images, in the current window
%     Or 4 full page images in figure windows 1-4
%     Or a single movie or slider
%     Or a single image in the axes specified by GRAPHHAND
%

%    Copyright (c) J. S. Marron 1999



%  Set Global variables
%
global GRAPHHAND A3DSLIMOV VSLIDERHAND VSTITLEHAND ;



%  Set defaults according to number of input arguments
%
if nargin == 1 ;    %  only 1 argument input, use default vitype
  ivitype = [1; 2; 1; 1; 0; 0] ;
else ;              %  xgrid was specified, use that
  ivitype = vitype ;

  if length(ivitype) == 1 ;   %  then need to add defaults
    ivitype = [ivitype; 2; 1; 1; 0; 0] ;
  elseif length(ivitype) == 2 ;   %  then need to add defaults
    if size(ivitype,1) == 1 ;   
      ivitype = [ivitype, 1, 1, 0, 0] ;  %  add to row vector
    elseif size(ivitype,2) == 1 ;   
      ivitype = [ivitype; 1; 1; 0; 0] ;  %  add to column vector
    end ;
  elseif length(ivitype) == 3 ;   %  then need to add default
    if size(ivitype,1) == 1 ;   
      ivitype = [ivitype, 1, 0, 0] ;  %  add to row vector
    elseif size(ivitype,2) == 1 ;   
      ivitype = [ivitype; 1; 0; 0] ;  %  add to column vector
    end ;
  elseif length(ivitype) == 4 ;   %  then need to add default
    if size(ivitype,1) == 1 ;   
      ivitype = [ivitype, 0, 0] ;  %  add to row vector
    elseif size(ivitype,2) == 1 ;   
      ivitype = [ivitype; 0; 0] ;  %  add to column vector
    end ;
  elseif length(ivitype) == 5 ;   %  then need to add default
    if size(ivitype,1) == 1 ;   
      ivitype = [ivitype, 0] ;  %  add to row vector
    elseif size(ivitype,2) == 1 ;   
      ivitype = [ivitype; 0] ;  %  add to column vector
    end ;
  end ;
end ;


if nargin <= 2 ;    %  at most 2 arguments input, use default var handling
  ivarp = [1,0] ;    %  Default
else ;
  ivarp = varp ;    %  value was specified, so use it
end ;


if nargin <= 3 ;    %  at most 3 arguments input, use default alpha
  ialpha = 0.05 ;    %  Default
else ;
  ialpha = alpha ;    %  value was specified, so use it
end ;


if nargin <= 4 ;    %  at most 4 arguments input, use default vhgp
  ivhgp = 0 ;
else ;              %  hgrid parameters were specified, use them

  if size(vhgp,2) > 1 ;    %  then have input row vector, so transpose
    ivhgp = vhgp' ;
  else ;
    ivhgp = vhgp ;
  end ;

end ;


if nargin <= 5 ;    %  at most 5 arguments input, use default vgp
  ivgp = 0 ;
else ;              %  grid parameters were specified, use them
  ivgp = vgp ;
end ;


if nargin <= 6 ;    %  at most 6 arguments input, use default bdryp
  ibdryp = 1 ;    %  Default
else ;
  ibdryp = bdryp ;    %  value was specified, so use it
end ;






%  Set parameters
%
arrowlength = 1.2 ;    %  in s2m1b, used 1.5
          %  length of arrows showing gradient direction
dotsize = 5 ;          %  from s2m7c, 5 was best for 64x64
          %  size of dots representing curvature
linelength = 0.2 ;     %  from s2m7c
          %  line length in pixel units  (0.8 gave "connected line segments")
          %  calibrated to be "1" for usual static set of four 64x64 images
zerosize = 2 ;         %  from s2m7c, used to use 3
          %  size of 0's representing low effective sample size

if ivitype(1) == 4 ;    %  then are doing streamlines, so
  ivitype(2) = 1 ;      %  reset to 1x1 single pixel analysis
end ;

if ivitype(2) == 2 ;
          %  if have chosen 2 x 2 blocks, and are not doing streamlines
          %  then need to modify size parameters for 2x2 blocks
  dotsize = 2 * dotsize ;
  zerosize = 2 * zerosize ;
          %  don't need similar change for arrowlength or linelength,
          %  since use one or the other
end ;


nrpt = 5 ;
          %  number of times to repeat movie
nrif = 4 ;
          %  number of times to Repeat Individual Movie Frames
          %    (controls speed of .mpeg version)
movopt = 1 ;
%  Movie options, currently using same as Cornean
%
%  This is MATLAB "help" to use for choosing parameters:
%	MPGWRITE(M, map, 'filename', options) Encodes M in MPEG format
% 	using the specified colormap and writes the result to the
% 	specified file.  The options argument is an optional vector of
% 	8 or fewer options where each value has the following meaning:
% 	1. REPEAT:
% 		An integer number of times to repeat the movie
% 		(default is 1).
% 	2. P-SEARCH ALGORITHM:
% 		0 = logarithmic	(fastest, default value)
% 		1 = subsample
% 		2 = exhaustive	(better, but slow)
% 	3. B-SEARCH ALGORITHM:
% 		0 = simple	(fastest)
% 		1 = cross2	(slightly slower, default value)
% 		2 = exhaustive	(very slow)
% 	4. REFERENCE FRAME:
% 		0 = original	(faster, default)
% 		1 = decoded	(slower, but results in better quality)
% 	5. RANGE IN PIXELS:
% 		An integer search radius.  Default is 10.
% 	6. I-FRAME Q-SCALE:
% 		An integer between 1 and 31.  Default is 8.
% 	7. P-FRAME Q-SCALE:
% 		An integer between 1 and 31.  Default is 10.
% 	8. B-FRAME Q-SCALE:
% 		An integer between 1 and 31.  Default is 25.
%movoutstr = '\matlab\steve\sss\sss1.mpg' ;
movoutstr = 'sss1.mpg' ;
          %  output string for movie mpeg file
          %  latest version writes to current directory



strmlntpb = 2.0 ; 
          %  for streamlines, average number of touches per pixel box
          %  This was chosen to be "right" for 16 x 16
          %  is appropriately rescaled below (after n,m are defined)
strmlnsf = .5 ;  
          %  for streamlines, step size factor, for multiplying pixel
          %  width to get step size
          %  This was chosen to be "right" for 16 x 16
          %  is appropriately rescaled below (after n,m are defined)
seed = 92374092 ;
rand('seed',seed) ;
          %  set Uniform random seed for use with streamlines





%  Set up colors for symbols (use scheme from Method 7 everywhere)
%
mcol = [1, 1, 0; ...
        1, .5, 0; ...
        1, 0, 0; ...
        .5, 0, 1; ...
        0, 0, .5] ; ...
%  yellow 
%  orange 
%  red 
%  purple 
%  dark blue 






%  detect whether image or density estimation data
%
if size(data,2) == 2 ;    %  Then is 2-d density estimation

  idatyp = 2 ;

  if ivitype(6) == 1 ;    %  then draw small ESS circles for density est.
    circleflag = 1 ;
  else ;    %  no circles for ESS (default)
    circleflag = 0 ;
  end ;

else ;   %  Then this is image data

  idatyp = 1 ;

  if ivitype(6) == 1 ;    %  then no small ESS circles for images
    circleflag = 0 ;
  else ;    %  draw circles for ESS (default)
    circleflag = 1 ;
  end ;

end ;




%  Set h grid stuff
%
if ivitype(3) == 1 ;    %  then are doing static plot

  if ivhgp == 0 ;   %  then use standard default h grid

    if idatyp == 1 ;    %  then are doing images
      ivhgp = [1; 8; 4] ;
    elseif idatyp == 2 ;    %  then are doing density est.
      ivhgp = [2; 16; 4] ;
    end ;

    iplottype = 1 ;
          %  1 - put 4x4 plot into current figure window
          %  2 - put 4 plots into figures 1,2,3,4
          %  3 - put single plot into current figure window
          %  4 - put single plot into input figure window
          %  5 - put single plot into input axes
  elseif length(ivhgp) == 1 ;   %  then use given bandwidth, and nh = 1
    ivhgp = [ivhgp; ivhgp; 1] ;

    if isempty(GRAPHHAND) ;
      iplottype = 3 ;
          %  1 - put 4x4 plot into current figure window
          %  2 - put 4 plots into figures 1,2,3,4
          %  3 - put single plot into current figure window
          %  4 - put single plot into input figure window
          %  5 - put single plot into input axes
    else ;
          %  put graphic in specified location

      if strcmp(get(GRAPHHAND(length(GRAPHHAND)),'Type'),'figure') ;
                                         %  then this is figure handle
        iplottype = 4 ;
          %  4 - put single plot into input figure window
      elseif strcmp(get(GRAPHHAND(length(GRAPHHAND)),'Type'),'axes') ;
                                         %  then this is axes handle
        iplottype = 5 ;
          %  5 - put single plot into input axes
      else ;
        iplottype = 3 ;
          %  3 - put single plot into current figure window
      end ;

    end ;


  elseif length(ivhgp) == 2 ;   %  then use given bandwidths, and nh = 4
    ivhgp = [ivhgp; 4] ;
    iplottype = 1 ;
          %  1 - put 4x4 plot into current figure window
  else ;

    if ivhgp(3) ~= 4 ;    %  then have bad nh
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      disp('!!!   Error from sss1.m:                  !!!') ;
      disp('!!!   For Static plot, must use nh = 4    !!!') ;
      disp('!!!   (or a single h)                     !!!') ;
      disp('!!!   Terminating sss1.m                  !!!') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

      return ;
          %  Stop all of sss1.m, and return to calling program

    end ;

    iplottype = 1 ;
          %  1 - put 4x4 plot into current figure window

  end ;

elseif  ivitype(3) == 2  |  ivitype(3) == 3  ;
                               %  then are doing dynamic or slider plot

  if ivhgp == 0 ;   %  then use standard default h grid

    if idatyp == 1 ;    %  then are doing images
      ivhgp = [1; 8; 25] ;
    elseif idatyp == 2 ;    %  then are doing density est.
      ivhgp = [2; 16; 25] ;
    end ;

  elseif length(ivhgp) == 2 ;   %  then use given bandwidths, and nh = 25
    ivhgp = [ivhgp; 25] ;
  elseif length(ivhgp) == 1 ;   %  single bandwidth?  give warning message

    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Warning from sss1.m:                !!!') ;
    disp('!!!   There is only a single bandwidth    !!!') ;
    disp('!!!   specified for this movie / slider   !!!') ;
    disp('!!!   Hit any key to Continue             !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

    pause ;

    ivhgp(2) = ivhgp(1) ;
    ivhgp(3) = 1 ;

  end ;

  if ~isempty(GRAPHHAND) ;    %  then have input handle
    if strcmp(get(GRAPHHAND(length(GRAPHHAND)),'Type'),'figure') ; 
                                        %  then this is figure handle
      iplottype = 4 ;
          %  1 - put 4x4 plot into current figure window
          %  2 - put 4 plots into figures 1,2,3,4
          %  3 - put single plot into current figure window
          %  4 - put single plot into input figure window
          %  5 - put single plot into input axes
    elseif strcmp(get(GRAPHHAND(length(GRAPHHAND)),'Type'),'axes') ;
                                         %  then this is axes handle

      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      disp('!!!   Warning from sss1.m:                 !!!') ;  
      disp('!!!   Cannot put this into axes,           !!!') ;
      disp('!!!   Will put in current figure instead   !!!') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

      iplottype = 3 ;
          %  3 - put single plot into current figure window
    else ;
      iplottype = 3 ;
          %  3 - put single plot into current figure window
    end ;

  else ;

    iplottype = 3 ;
          %  3 - put single plot into current figure window

  end ;

end ;

hmin = ivhgp(1) ;
hmax = ivhgp(2) ;
nh = ivhgp(3) ;
if nh == 1 ;
  vh = hmin ;
else ;
  vh = logspace(log10(hmin),log10(hmax),nh) ;
end ;

if ivitype(3) == 2 ;    %  then are making movie

  vmorder = [(1:nh),((nh-1):-1:2)] ;
  vmorder = vec2mat(vmorder,nrif) ;
  vmorder = reshape(vmorder,1,size(vmorder,2)*nrif) ;
          %  replicate frame number, NRIF times

elseif ivitype(3) == 3 ;    %  then are doing a slider plot

  %  first set up global variables
  %
  global SLIMOV VLOG10H CURRENTH SLIDERHAND STITLEHAND SSSSTR ;
          %  SLIMOV - slider's movie matrix
          %  VLOG10H - vector of log10 bandwidths
          %  CURRENTH - current value of bandwidth for slider
          %  SLIDERHAND - Handle for Slider (to get value)
          %  STITLEHAND - Handle for Title of Slider (showing bandwidth value)
          %  SSSSTR - String for type of SSS that is running

  if iplottype == 4 ;    %  then use input figure handle
    sss1slider(vh,ivitype) ;
          %  Set up GUI window
  else ;    %  then use current figure
    clf ;
          % clear current figure window
    sss1slider(vh,ivitype) ;
          %  Set up GUI window
  end ;

end ;




%  Do preliminary calculations
%
if idatyp == 2 ;    %  Then is 2-d density estimation

  %  Set x grid stuff
  %
  if length(ivgp) == 1 ;   %  then use standard default x grid
    xmin = min(data(:,1)) ;
    xmax = max(data(:,1)) ;
    xrange = xmax - xmin ;
    xmin = xmin - 0.05 * xrange ;
    xmax = xmax + 0.05 * xrange ;
    nxg = 64 ;
    ymin = min(data(:,2)) ;
    ymax = max(data(:,2)) ;
    yrange = ymax - ymin ;
    ymin = ymin - 0.05 * yrange ;
    ymax = ymax + 0.05 * yrange ;
    nyg = 64 ;

    if ivgp == 1 ;     %  then use nearest neighbor binning
      bintype = 1 ;    
    else ;
      bintype = 2 ;    
    end ;

  else ;
    xmin = ivgp(1) ;
    xmax = ivgp(2) ;
    nxg = ivgp(3) ;
    ymin = ivgp(4) ;
    ymax = ivgp(5) ;
    nyg = ivgp(6) ;

    if (round(nxg) ~= nxg)  |  (round(nyg) ~= nyg) ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      disp('!!!   Warning from SSS1:  input grid sizes        !!!') ;
      disp('!!!   are not integers, will use rounded values   !!!') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

      nxg = round(nxg)
      nyg = round(nyg)
    end ;

    if length(ivgp) == 6 ;     %  then use linear binning
      bintype = 2 ;    
    else ;
      bintype = 1 ;    
    end ;

  end ;



  %  Bin the data
  %
  data = sss1bin(data,xmin,xmax,nxg,ymin,ymax,nyg,bintype,ibdryp) ;
          %  Note data are converted from Cartesian to 
          %  Matrix coordinates inside this
          


else ;                    %  Then assume image data


  if ibdryp == 1 ;    %  then do boundary adjustment, based on mean
    data = data - mean(mean(data)) ;
  elseif ibdryp == 2 ;    %  then do boundary adjustment, based on median
    data = data - median(median(data)) ;
  end ;


end ;





%  Set additional parameters
%
m = size(data,2) ;
          %  number of cols of data matrix 
          %  (i.e. j-values in matrix coordinate system)
n = size(data,1) ;
          %  for images: number of rows of data matrix 
          %  (i.e. i-values in matrix coordinate system)
maxnm = max([n; m]) ;
dotsize = (64 / maxnm) * dotsize ;
zerosize = (64 / maxnm) * zerosize ;
          %  resize, since original calibration was for n = 64,
          %  i.e. put dot size on pixel scale
          %  Note:  need to do this with dots, and not with arrows,
          %      since, arrows are naturally (by their construction)
          %      scaled to pixel size, but dots are not.

strmlntpb = strmlntpb * (16 / maxnm) ; 
          %  for streamlines, average number of touches per pixel box
strmlnsf = .5 * (maxnm / 16) ;  
          %  for streamlines, step size factor, for multiplying pixel
          %  width to get step size




%  Check for reasonable image sizes
%
if (ivitype(5) == 0) & (nh == 4) ;   %  then do sample size warnings and prompts
                                     %  note: don't do this if already
                                     %  have single page

  if ivitype(1) == 1  |  ivitype(1) == 3 ;    %  then have some arrows

    if ((ivitype(2) == 1) & (maxnm > 64)) |...
            ((ivitype(2) == 2) & (maxnm > 128)) ;
                           %  if working with 1x1 single pixels and n > 64
                           %  or working on 2x2 blocks and n > 128
                           %  then arrows won't look good, so give warning

      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      disp('!!!   Warning from sss1.m:                        !!!!') ;
      disp(['!!!   Image Size = ' num2str(maxnm)]) ;
      disp('!!!   is too large for Good Resolution of Arrows   !!!') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      innum = input(['!!!   Choose from: \n' ...
             '!!!      1.  Continue anyway (press Enter)\n' ...
             '!!!      2.  Make full page plots (press 2 + Enter)\n' ...
             '!!!      3.  Quit processing  (press 3 + Enter)\n' ...
             '!!!            (e.g. go to movie, or single h)\n']) ;

      if innum == 2 ;

        iplottype = 2 ;  
          %  1 - put 4x4 plot into current figure window
          %  2 - put 4 plots into figures 1,2,3,4
          %  3 - put single plot into current figure window
          %  4 - put single plot into input figure window
          %  5 - put single plot into input axes

      elseif innum == 3 ;    %  Then want to quit executing

        return ;
          %  Stop all of sss1.m, and return to calling program
      end ;

    end ;

  end ;

elseif ivitype(5) == 2 ;    %  then go to single page when needed

  if (ivitype(1) == 1  |  ivitype(1) == 3) & ... 
     (((ivitype(2) == 1) & (maxnm > 64)) |...
            ((ivitype(2) == 2) & (maxnm > 128))) ;
                           %  if have some arrows and
                           %  are working with 1x1 single pixels and n > 64
                           %  or working on 2x2 blocks and n > 128
                           %  then arrows won't look good, so go to full page

    iplottype = 2 ;  
          %  1 - put 4x4 plot into current figure window
          %  2 - put 4 plots into figures 1,2,3,4
          %  3 - put single plot into current figure window
          %  4 - put single plot into input figure window
          %  5 - put single plot into input axes

  end ;

end ;



if iplottype == 1 ;
  clf ;
          %  clear current figure window
end ;




%  Loop through bandwidths, i.e. pictures
%
for ih = 1:nh ;

  h = vh(ih) ;


  if ivitype(4) == 1 ;    %  then display progress messages
    disp(['    SSS1:  Working on h = ' num2str(h)]) ;
  end ;



  %  First do smooths and Effective Sample Size
  %
  if idatyp == 1 ;    %  then are doing regression

    if ivitype(1) == 0 ;   %  then are doing only smooths

      fh = sssfh(data,h,-1) ;

    elseif ivarp(1) == 0 ;   %  then use input sigma^2 (noise level of data)

      [fh,ess] = sssfh(data,h) ;

      sig2 = ivarp(2) ;
      sig2 = sig2 * ones(n,m) ;

    elseif ivarp(1) == 1 ;    %  then work with estimated variance

      [fh,ess,sig2] = sssfh(data,h,ivarp(2)) ;

    end ;

  elseif idatyp == 2 ;    %  then are doing density estimation

    if ivitype(1) == 0 ;   %  then are doing only smooths

      fh = sssfhd(data,h,-1) ;

    else ;    %  then work with estimated variance

      [fh,ess,sig2] = sssfhd(data,h) ;

    end ;

  end;

  if ivitype(4) == 1 ;    %  then display progress messages
    disp('    SSS1:    fh, ESS and sig2 finished') ;
  end ;




  %  Plot current smooth as gray level image
  %
  if iplottype == 1 ;       %  1 - put 4x4 plot into current figure window
    subplot(2,2,ih) ;
  elseif iplottype == 2 ;   %  2 - put 4 plots into figures 1,2,3,4
    figure(ih) ;
    clf ;
    subplot(1,1,1) ;
  else ;                    %  3 - put single plot into current figure window
                            %  4 - put single plot into input figure window
                            %  5 - put single plot into input axes

    if ~(ivitype(3) == 3) ;     %  the don't have a slider plot
      subplot(1,1,1) ;
    end ;

  end ;


  colormap(gray) ;
    mfh = fh - min(min(fh)) ;
    mfh = mfh / max(max(mfh)) ;
    mfh = mfh * size(gray,1) ;
          %  version of fh, mapped to [0, #gray levels]
    image([1,m],[1,n],mfh) ;
      if ivitype(3) == 1 ;    %  then are doing static plot
        title(['h = ' num2str(h)]) ;
      end ;


  if idatyp == 1 ;    %   then are doing regression,
                    %   so want each pixel to be a square
    axis('image') ;

  elseif idatyp == 2 ;    %  then are doing density estimation, 
                          %  so relabel axes, and
                          %  want to set aspect ratio for density
                          %  with respect to 2-d Lebesgue measure

    axis('image') ;

    imxticks = get(gca,'XTick') ;
    imyticks = get(gca,'YTick') ;
          %  Value for tick marks, in image coordinates

    caxticks = (xmax - xmin) * (imxticks / nxg) + xmin ;
    cayticks = (ymin - ymax) * (imyticks / nyg) + ymax ;
          %  Value for tick marks, in cartesian coordinates
          %  (recall y axis is reversed)

    caxticklabels = num2str(caxticks',3) ;
    cayticklabels = num2str(cayticks',3) ;
          %  Convert these to strings for use in graphics

    set(gca,'XTickLabel',caxticklabels) ;
    set(gca,'YTickLabel',cayticklabels) ;

  end ;



  if ivitype(3) == 3 ;    %  then are doing a slider plot,
                              %      so update text and slider position

    set(SLIDERHAND,'Value',log10(h)) ;

    CURRENTH = h ;
      tstr = [SSSSTR ', h = ' num2str(CURRENTH)] ;
    set(STITLEHAND,'String',tstr) ;

  end ;


  if ivitype(1) ~= 0 ;   %  then want to add some SSS parts

    if circleflag == 1 ;    %  then want to draw small ESS circles

      %  Do Small Effective Sample Size Calculations
      %
      smallessflag = ess < 5 ;
          %  ones where ess < 5 ;
      if ivitype(2) == 1 ;    %  then are working with single pixels
        nsmalless = sum(sum(smallessflag)) ;
          %  number of pixels with small ess
        allsmalless = (nsmalless == n * m) ;
          %  1 when every pixel has small ess
        nosmalless = (nsmalless == 0) ;
          %  1 when no pixel has small ess
        smallessflagdec = [] ;
          %  can pass this into sss1ssp.m
      elseif ivitype(2) == 2 ;    %  then are working with 2x2 blocks
        mo2 = floor(m/2) ;
        no2 = floor(n/2) ;
        smallessflagdec = conv2(smallessflag,ones(2,2),'valid') ;
          %  do 2x2 simple sum
        smallessflagdec = smallessflagdec((1:2:(2*no2)),(1:2:(2*mo2))) ;
          %  reduce to just block centers
        smallessflagdec = smallessflagdec > 0 ;
          %  puts a one if any cell has small ess
        nsmalless = sum(sum(smallessflagdec)) ;
          %  number of 2x2 blocks with small ess
        allsmalless = (nsmalless == no2 * mo2) ;
          %  1 when every 2x2 block has small ess
        nosmalless = (nsmalless == 0) ;
          %  1 when no 2x2 block has small ess
      end ;


    else ;    %  Don't draw any smalless circles

      nsmalless = 0 ;
          %  number of pixels with small ess
      allsmalless = logical(0) ;
          %  1 when every pixel has small ess
      nosmalless = logical(1) ;
          %  1 when no pixel has small ess
      smallessflag = logical(zeros(n,m)) ;

      if ivitype(2) == 1 ;    %  then are working with single pixels
        smallessflagdec = [] ;
          %  can pass this into sss1ssp.m
      elseif ivitype(2) == 2 ;    %  then are working with 2x2 blocks
        mo2 = floor(m/2) ;
        no2 = floor(n/2) ;
        smallessflagdec = logical(zeros(no2,mo2)) ;
          %  puts a one if any cell has small ess
      end ;


    end ;




    if allsmalless ;          %  then Effective Sample Size is nowhere
                              %  large enough, so don't do hard calculations
                              %  but overlay 0's everywhere


      sss1asp(n,m,ivitype(2),zerosize) ;



    else ;       %  then go ahead with derivative calculations


      if ivitype(1) == 1  |  ivitype(1) == 3  |  ivitype(1) == 4 ; 
                                      %  Then do gradient type calcs


        [vsigg, marrow] = sssgrad(data,h,ivitype(2),ess,sig2,ialpha,idatyp) ;
          %  vitype(2) is "pixel parameter" (1 or 2)


        %  eliminate arrows in small sample size areas
        %
        if ivitype(2) == 1 ;    %  then are working with single pixels
          vsigg = vsigg & ~reshape(smallessflag,n*m,1) ;
            %  ones where have arrows & enough data
        elseif ivitype(2) == 2 ;    %  then are working with 2x2 blocks
          vsigg = vsigg & ~reshape(smallessflagdec,no2*mo2,1) ;
            %  ones where have arrows & enough data
        end ;
        flag = (vsigg(marrow(:,1)) == 1) ;
          %  ones in rows of marrow, that will keep
        marrow = marrow(flag,:) ;


        if ivitype(4) == 1 ;    %  then display progress messages
          disp('    SSS1:    gradient calculations finished') ;
        end ;

      else ;

        vsigg = [] ;
            %  create this for passing to sss1dot.m

      end ;




      if ivitype(1) == 2  |  ivitype(1) == 3 ;  %  Then do curvature type calcs

        [vsigc, mdot] = ssscurv(data,h,vitype(2),ess,sig2,ialpha,idatyp) ;
          %  vitype(2) is "pixel parameter" (1 or 2)


        %  eliminate dots in small sample size areas
        %
        if ivitype(2) == 1 ;    %  then are working with single pixels
          vsigc = vsigc & ~reshape(smallessflag,n*m,1) ;
            %  ones where have dots & enough data
        elseif ivitype(2) == 2 ;    %  then are working with 2x2 blocks
          vsigc = vsigc & ~reshape(smallessflagdec,no2*mo2,1) ;
            %  ones where have dots & enough data
        end ;
        flag = (vsigc(mdot(:,1)) == 1) ;
          %  ones in rows of mdot, that will keep
        mdot = mdot(flag,:) ;


        if ivitype(4) == 1 ;    %  then display progress messages
          disp('    SSS1:    curvature calculations finished') ;
        end ;

      end ;




      %  First set up colors (both arrows and dots)
      %
      if ivitype(2) == 1 ;    %  then are working with single pixels
        colormask = zeros(n,m) ;
      elseif ivitype(2) == 2 ;    %  then are working with 2x2 blocks
        mo2 = floor(m/2) ;
        no2 = floor(n/2) ;
        colormask = zeros(no2,mo2) ;
      end ;
          %  initially set to "no significant curvature"

      if ivitype(1) == 2  |  ivitype(1) == 3 ;  
                    %  Then update colormap with curvature colors

        ndot = sum(vsigc) ;    
          %  number of dots to draw

        if ndot > 0 ;    %  then have some dots

          colormask(mdot(:,1)) = mdot(:,2) ;
            %  replace entry numbers in mdot(:,1), 
            %  with sig codes in mdot(:,2)

        end ;

      end ;




      %  Draw arrows
      %
      if ivitype(1) == 1  |  ivitype(1) == 3 ;  
                    %  Then have done gradient calculations, so plot arrows

        sss1arrp(marrow,vsigg,colormask,mcol,n,m,ivitype(2), ...
                                              linelength, arrowlength) ;

      end ;




      %  Draw dots
      %
      if ivitype(1) == 2  |  ivitype(1) == 3 ;  
                    %  Then have done curvature calculations

        sss1dotp(mdot,vsigc,vsigg,colormask,mcol,n,m,ivitype(1),...
                                                     ivitype(2),dotsize) ;
      end ;




      %  Draw Streamlines
      %
      if ivitype(1) == 4 ;  
                    %  Then have done curvature calculations

        sss1strmln(marrow,vsigg,n,m,strmlntpb,strmlnsf) ;

      end ;





      %  Put down Small ESS Zeros as needed
      %
      if nosmalless == 0 ;     %  then have some small ess parts to show

        sss1ssp(smallessflag,smallessflagdec,n,m,ivitype(2),zerosize) ;

      end ;

    end ;


  end ;




  if  ivitype(3) == 2  |  ivitype(3) == 3  ;
                               %  then are doing dynamic or slider plot

    %  Store current image as a frame
    %
    sssmov(:,ih) = getframe(gca) ;


  end ;



end ;






%  Play movies or construct slider  (if needed)
%
if ivitype(3) == 2 ;    %  then make movie

  sssmov = sssmov(:,vmorder') ;
          %  reorder frames, to loop back to beginning

  movie(sssmov,nrpt) ;
          %  Play movie on screen

  mpgwrite(sssmov,gray,movoutstr,movopt) ;
          %  save as mpeg file



elseif ivitype(3) == 3  ;   % then make slider plot

  SLIMOV = sssmov ;
          %  make these global for use by slider

  if iplottype == 4 ;
    A3DSLIMOV = cat(3,A3DSLIMOV,SLIMOV) ;
    VSLIDERHAND = [VSLIDERHAND, SLIDERHAND] ;
    VSTITLEHAND = [VSTITLEHAND, STITLEHAND] ;
  end ;

end ;



