PURPOSE:   (one line only)
  Optimal image subtraction

  STAMP - Definition - This term refers to a localized section of the array
    that contains an image of point-source.  This could either be an
    actual sub-array or a description of a sub-array in an larger image.

  CCD data processing
  image - Image to be processed.  If this is a string it is taken to be
             a file name that will be read.  The other option is to provide
             the array to be processed.
  reference - Reference image to be subtracted from image.  Just as with
                 image, the input can be a string or an array.
  CONSTPHOT - Flag, if set turns on the constant photometric ratio constraint
  CHISQSTAT - Threshold for excluding stamps that don't fit well.
                 Default=0.0 (ie., no exclusions)
  GAUSSIAN  - Keyword that controls the basis set used.  If not provided,
                a Delta-function basis set is used.  To use a gaussian
                function basis set either give the name of the basis set
                to use or privide the description yourself.  If you provide
                'astier' as the name, you will get the following:
              gaussian={nc:3, degree:[6,4,2], sigmas:[0.7,1.5,2.0]}
              You can provide your own description by providing a
              similar structure.  Tag nc is the number of components.
              The length of degree and sigmas must match nc.
  PATH - Directory location to find the input images.  Default is the
           current directory.
  IMGHDR - String array containing FITS header information for input image.
             This is used only when the input image is provided directly,
             rather than to be read from a file.  This array is used as
             a starting point for the saved output image (if desired).
  FWHM - FWHM of the input image.  Used only if image provided
             directly and the header is not.  This value is
             otherwise read from the header from the 'SEEING'
             keyword in the header.
  MAXPHOTSIG- Maximum DN value for a useful signal.  Any source with a peak
                above this level is passed over.  Default=60000.0 DN
  MINFLUX - Threshold for sources to be included in the reference PSFs
               to support the image subtraction.  This value is relative
               to the brightest source found in the image.  Default=0.25
  DEGREE - Degree of polynomial used to fit for the space varying kernel.
             Default=0 (constant kernel)
  DELP   - Grid size in pixels for recalculating the space-varying kernel.
             Default=1 (doesn't matter for a constant kernel)
  OUTPATH - String, name of directory to save data in.
                Default = current directory.
                NOTE: PATH and OUTPATH must not be the same.  If they are
                  then the SAVEIT flag will be ignored and the output will
                  not be generated.
  SILENT - Flag that will determine whether various relevant information will
           be printed. Default=0. By default, it will print the statistics
  NODISPLAY - Flag that will determine whether various relevant plots will be
              created. Default=0. By default, the plots will be created.
  SEED    - value for the seed used by the random call. Default is undefined
              causing random to use the system time as the seed
  diffimage - the subtracted image 
  CREFERENCE - the convolved reference image
  FAILED - boolean indicating whether or not the subtraction failed
  DIFFINFO - anonymous structure containing lots of useful information about
                the difference image and what happened along the way.
                This structure is undefined if the output keyword FAILED is set
  Note that the input to this routine is expected to already be registered
    and be of the same size.  The routine, dewarp, is usually used to
    do this registration step.

  This routine struggles as the image FWHM gets a lot larger than the
    template.  The execution time goes up significantly as well.

   Internal helper routines:
   Written by SwRI clinic team, 2009/11/09, optimized version derived from
     the original code written by Patrick Miller. (ref)
   2010/03/05, FS, added IMGSKY and REFSKY keywords
   2010/03/12, FS, added SILENT, NODISPLAY, and SEED keywords
   2010/04/21, FS, added CREFERENCE, FAILED keywords
   2010/04/28, FS, changed shift calls in ois_coef to sshift2d with edge_zero
   2010/06/10, MWB, fixed a display bug with the number of stamps changing.
                     Also changed default on DEGREE to 0.
   2010/09/23, MWB, maxphotsig was not getting all the way to the findsrc
                     call that does the stamp selection.  This has been
                     changed.  Also, kernel size determination was tweaked
   2010/10/13, MWB, rework of the image/reference flux ratio calculation
   2010/10/21, MWB, fixed a bug where the FWHM value read from the header
                      was not used properly.
   2010/10/22, MWB, cleaded up ois_stamps routine.  Flux ratio and fwhm
                      determinations are much better now.
   2010/11/02, MWB, added DIFFINFO output keyword
   2010/11/03, MWB, minor modification to ois_convref to hold the output
                      image type to float (was returning a double).  This
                      ensures the final difference image is a float.
   2011/04/12, MWB, fixed a bug in stamp selection for a rare case when
                      the brightest good stamp is too close to the edge
                      of the image.
   2011/12/28, MWB, Changed usage of MAXPHOTSIG.  Now, this level is used
                      just during the stamp selection process.  No change
                      is made during the actual differencing process.
   2012/01/24, MWB, Changed FWHM filtering on stamps.  It was a 2-sigma cut
                      and this has been relaxed to a 3-sigma cut.  Also,
                      one case that generated a warning and tried to
                      run the difference anyway was changed to return with
                      a failure.
   2012/03/10, MWB, Found and fixed a minor bug in ois_fluxcheck.  It would
                      sometimes cause the differencing to fail if the first
                      stamp had already been marked bad.  This routine was
                      changed to use the first good stamp.
  2014/02/11, MWB, added more data type options for SEED
  2016/05/02, MWB, integrated new ois_coeff routine from Steve Hartung that
                      fixes a bug related to high-order PSF variations