Initial commit

21 files changed, 6040 insertions, 0 deletions

diff --git a/pkg/images/imfit/Revisions b/pkg/images/imfit/Revisions
new file mode 100644
index 00000000..7ff95abd
--- /dev/null
+++ b/pkg/images/imfit/Revisions
@@ -0,0 +1,2025 @@
+.help revisions Jan97 images.imfit
+.nf
+===============================
+Package Reorganization
+===============================
+
+pkg/images/imarith/t_imsum.x
+pkg/images/imarith/t_imcombine.x
+pkg/images/doc/imsum.hlp
+pkg/images/doc/imcombine.hlp
+    Provided options for USHORT data.  (12/10/96, Valdes)
+
+pkg/images/imarith/icsetout.x
+pkg/images/doc/imcombine.hlp
+    A new option for computing offsets from the image WCS has been added.
+    (11/30/96, Valdes)
+
+pkg/images/imarith/t_imcombine.x
+pkg/images/imarith/icombine.gx
+    Changed the error checking to catch additional errors relating to too
+    many files.  (11/12/96, Valdes)
+
+pkg/images/imarith/icsort.gx
+    There was an error in the ic_2sort routine when there are exactly
+    three images that one of the explicit cases did not properly keep
+    the image identifications.  See buglog 344.  (8/1/96, Valdes)
+
+pkg/images/filters/median.x
+    The routine mde_yefilter was being called with the wrong number of
+    arguments.
+    (7/18/96, Davis)
+
+pkg/images/imarith/t_imcombine.x
+pkg/images/imarith/icombine.gx
+pkg/images/imarith/icimstack.x +
+pkg/images/imarith/iclog.x
+pkg/images/imarith/mkpkg
+pkg/images/doc/imcombine.hlp
+    The limit on the maximum number of images that can be combined, set by
+    the maximum number of logical file descriptors, has been removed.  If
+    the condition of too many files is detected the task now automatically
+    stacks all the images in a temporary image and then combines them with
+    the project option.
+    (5/14/96, Valdes)
+
+pkg/images/geometry/xregister/rgxfit.x
+    Changed several Memr[] references to Memi[] in the rg_fit routine.
+    This bug was causing a floating point error in the xregister task
+    on the Dec Alpha if the coords file was defined, and could potentially
+    cause problems on other machines.
+    (Davis, April 3, 1996)
+
+pkg/images/geometry/t_geotran.x
+pkg/images/geometry/geograph.x
+pkg/images/doc/geomap.hlp
+     Corrected the definition of skew in the routines which compute a geometric
+     interpretation of the 6-coefficient fit, which compute the coefficients
+     from the geometric parameters, and in the relevant help pages.
+     (2/19/96, Davis)
+
+pkg/images/median.par
+pkg/images/rmedian.par
+pkg/images/mode.par
+pkg/images/rmode.par
+pkg/images/fmedian.par
+pkg/images/frmedian.par
+pkg/images/fmode.par
+pkg/images/frmode.par
+pkg/images/doc/median.hlp
+pkg/images/doc/rmedian.hlp
+pkg/images/doc/mode.hlp
+pkg/images/doc/rmode.hlp
+pkg/images/doc/fmedian.hlp
+pkg/images/doc/frmedian.hlp
+pkg/images/doc/fmode.hlp
+pkg/images/doc/frmode.hlp
+pkg/images/filters/t_median.x
+pkg/images/filters/t_rmedian.x
+pkg/images/filters/t_mode.x
+pkg/images/filters/t_rmode.x
+pkg/images/filters/t_fmedian.x
+pkg/images/filters/t_frmedian.x
+pkg/images/filters/t_fmode.x
+pkg/images/filters/t_frmode.x
+    Added a verbose parameter to the median, rmedian, mode, rmode, fmedian,
+    frmedian, fmode, and frmode tasks. (11/27/95, Davis)
+
+pkg/images/geometry/doc/geotran.hlp
+    Fixed an error in the help page for geotran. The default values for
+    the xscale and yscale parameters were incorrectly listed as INDEF,
+    INDEF instead of 1.0, 1.0. (11/14/95, Davis)
+
+pkg/images/imarith/icpclip.gx
+    Fixed a bug where a variable was improperly used for two different
+    purposes causing the algorithm to fail (bug 316).  (10/19/95, Valdes)
+
+pkg/images/doc/imcombine.hlp
+    Clarified a point about how the sigma is calculated with the SIGCLIP
+    option.  (10/11/95, Valdes)
+
+pkg/images/imarith/icombine.gx
+    To deal with the case of readnoise=0. and image data which has points with
+    negative mean or median and very small minimum readnoise is set
+    internally to avoid computing a zero sigma and dividing by it.  This
+    applies to the noise model rejection options.  (8/11/95, Valdes)
+
+pkg/images/frmedian.hlp
+pkg/images/frmode.hlp
+pkg/images/rmedian.hlp
+pkg/images/rmode.hlp
+pkg/images/frmedian.par
+pkg/images/frmode.par
+pkg/images/rmedian.par
+pkg/images/rmode.par
+pkg/images/filters/frmedian.h
+pkg/images/filters/frmode.h
+pkg/images/filters/rmedian.h
+pkg/images/filters/rmode.h
+pkg/images/filters/t_frmedian.x
+pkg/images/filters/t_frmode.x
+pkg/images/filters/t_rmedian.x
+pkg/images/filters/t_rmode.x
+pkg/images/filters/frmedian.x
+pkg/images/filters/frmode.x
+pkg/images/filters/rmedian.x
+pkg/images/filters/rmode.x
+pkg/images/filters/med_utils.x
+    Added new ring median and modal filtering tasks frmedian, rmedian,
+    frmode, and rmode to the images package.
+    (6/20/95, Davis)
+
+pkg/images/fmedian.hlp
+pkg/images/fmode.hlp
+pkg/images/median.hlp
+pkg/images/mode.hlp
+pkg/images/fmedian.par
+pkg/images/fmode.par
+pkg/images/median.par
+pkg/images/mode.par
+pkg/images/filters/fmedian.h
+pkg/images/filters/fmode.h
+pkg/images/filters/median.h
+pkg/images/filters/mode.h
+pkg/images/filters/t_fmedian.x
+pkg/images/filters/t_fmode.x
+pkg/images/filters/t_median.x
+pkg/images/filters/t_mode.x
+pkg/images/filters/fmedian.x
+pkg/images/filters/fmode.x
+pkg/images/filters/median.x
+pkg/images/filters/mode.x
+pkg/images/filters/fmd_buf.x
+pkg/images/filters/fmd_hist.x
+pkg/images/filters/fmd_maxmin.x
+pkg/images/filters/med_buf.x
+pkg/images/filters/med_sort.x
+    Added minimum and maximum good data parameters to the fmedian, fmode,
+    median, and mode filtering tasks.  Removed the 64X64 kernel size limit
+    in the median and mode tasks.  Replaced the common blocks with structures
+    and .h files.
+    (6/20/95, Davis)
+
+pkg/images/geometry/t_geotran.x
+pkg/images/geometry/geotran.x
+pkg/images/geometry/geotimtran.x
+    Fixed a bug in the buffering of the x and y coordinate surface interpolants
+    which can cause a memory corruption error if, nthe nxsample or nysample
+    parameters are > 1, and the nxblock or nyblock parameters are less
+    than the x and y dimensions of the input image. Took the opportunity
+    to clean up the code.
+    (6/13/95, Davis)
+
+=======
+V2.10.4
+=======
+
+pkg/images/geometry/t_geomap.x
+    Corrected a harmless typo in the code which determines the minimum
+    and maximum x values and improved the precision of the test when the
+    input is double precision.
+    (4/18/95, Davis)
+
+pkg/images/doc/fit1d.hlp
+    Added a description of the interactive parameter to the fit1d help page.
+    (4/17/95, Davis)
+
+pkg/images/imarith/t_imcombine.x
+    If an error occurs while opening an input image header the error
+    recovery will close all open images and then propagate the error.
+    For the case of running out of file descriptors with STF format
+    images this will allow the error message to be printed rather
+    than the error code.  (4/3/95, Valdes)
+
+pkg/images/geometry/xregister/t_xregister.x
+    Added a test on the status code returned from the fitting routine so
+    the xregister tasks does not go ahead and write an output image when
+    the user quits the task in in interactive mode.
+    (3/31/95, Davis)
+
+pkg/images/imarith/icscale.x
+pkg/images/doc/imcombine.hlp
+    The behavior of the weights when using both multiplicative and zero
+    point scaling was incorrect; the zero levels have to account for
+    the scaling.  (3/27/95, Valdes)
+
+pkg/images/geometry/xregister/rgxtools.x
+    Changed some amovr and amovi calls to amovkr and amovki calls.
+    (3/15/95, Davis)
+
+pkg/images/geometry/t_imshift.x
+pkg/images/geometry/t_magnify.x
+pkg/images/geometry/geotran.x
+pkg/images/geometry/xregister/rgximshift.x
+    The buffering margins set for the bicubic spline interpolants were
+    increased to improve the flux conservation properties of the interpolant
+    in cases where the data is undersampled. (12/6/94, Davis)
+
+pkg/images/xregister/rgxbckgrd.x
+    In several places the construct array[1++nx-wborder] was being used
+    instead of array[1+nx-wborder]. Apparently caused by a typo which
+    propagated through the code, the Sun compilers did not catch this, but
+    the IBM/RISC6000 compilers did. (11/16/94, Davis)
+
+
+pkg/images/xregister.par
+pkg/images/doc/xregister.hlp
+pkg/images/geometry/xregister/t_xregister.x
+pkg/images/geometry/xregister/rgxcorr.x
+pkg/images/geometry/xregister/rgxicorr.x
+pkg/images/geometry/xregister/rgxcolon.x
+pkg/images/geometry/xregister/rgxdbio.x
+    The xregister task was modified to to write the output shifts file
+    in either text database format (the current default)  or in simple text
+    format. The change was made so that the output of xregister could
+    both be edited more easily by the user and be used directly with the
+    imshift task. (11/11/94, Davis)
+
+pkg/images/imfit/fit1d.x
+    A Memc in the ratio output option was incorrectly used instead of Memr
+    when the bug fix of 11/16/93 was made.  (10/14/94, Valdes)
+
+pkg/images/geometry/xregister/rgxcorr.x
+    The procedure rg_xlaplace was being incorrectly declared as an integer
+    procedure.
+    (8/1/94, Davis)
+
+pkg/images/geometry/xregister/rgxregions.x
+    The routine strncmp was being called (with a missing number of characters
+    argument) instead of strcmp. This was causing a bus error under solaris
+    but not sun os whenever the user set regions to "grid ...". (7/27/94 LED)
+
+pkg/images/tv/imexaine/ierimexam.x
+    The Gaussian fitting can return a negative sigma**2 which would cause
+    an FPE when the square root is taken.  This will only occur when
+    there is no reasonable signal.  The results of the gaussian fitting
+    are now set to INDEF if this unphysical result occurs.  (7/7/94, Valdes)
+
+pkg/images/geometry/geofit.x
+    A routine expecting two char arrays was being passed two real arrays
+    instead resulting in a segmentation violation if calctype=real
+    and reject > 0.
+    (6/21/94, Davis)
+
+pkg/images/imarith/t_imarith.x
+    IMARITH now deletes the CCDMEAN keyword if present.  (6/21/94, Valdes)
+
+pkg/images/imarith/icaclip.gx
+pkg/images/imarith/iccclip.gx
+pkg/images/imarith/icpclip.gx
+pkg/images/imarith/icsclip.gx
+    1.	The restoration of deleted pixels to satisfy the nkeep parameter
+	was being done inside the iteration loop causing the possiblity
+	of a non-terminating loop; i.e. pixels are rejected, they are
+	restored, and the number left then does not statisfy the termination
+	condition.  The restoration step was moved following the iterative
+	rejection.
+    2.	The restoration was also incorrectly when mclip=no and could
+	lead to a segmentation violation.
+    (6/13/94, Valdes)
+
+pkg/images/geometry/xregister/rgxicorr.x
+    The path names to the xregister task interactive help files was incorrect.
+    (6/13/94, Davis)
+
+pkg/images/imarith/iccclip.gx
+pkg/images/imarith/icsclip.gx
+    Found and fixed another typo bug.  (6/7/94, Valdes/Zhang)
+
+pkg/images/imarith/icscale.x
+    The sigma scaling flag, doscale1, would not be set in the case of
+    a mean offset of zero though the scale factors could be different.
+    (5/25/94, Valdes/Zhang)
+
+pkg/images/imarith/icsclip.gx
+    There was a missing line: l = Memi[mp1].  (5/25/94, Valdes/Zhang)
+
+pkg/images/imarith/icaclip.gx
+pkg/images/imarith/iccclip.gx
+pkg/images/imarith/icpclip.gx
+pkg/images/imarith/icsclip.gx
+    The reordering step when a central median is used during rejection
+    but the final combining is average was incorrect if the number
+    of rejected low pixels was greater than the number of pixel
+    number of pixels not rejected.  (5/25/94, Valdes)
+
+pkg/images/geometry/t_geotran.x
+    In cases where there was no input geomap database, geotran was
+    unnecessarily  overiding the size of the input image requested by the
+    user if the size of the image was bigger than the default output size
+    (the size of the output image which would include all the input image
+    pixels is no user shifts were applied).
+    (5/10/94, Davis)
+
+pkg/images/imarith/icscale.x
+pkg/images/imarith/t_imcombine.x
+    1.  There is now a warning error if the scale, zero, or weight type
+	is unknown.
+    2.  An sfree was being called before the allocated memory was finished
+	being used.
+    (5/2/94, Valdes)
+
+pkg/images/tv/imexaine/ierimexam.x
+    For some objects the moment analysis could fail producing a floating
+    overflow error in imexamine, because the code was trying to use
+    INDEF as the initial value of the object fwhm.  Changed the gaussian
+    fitting code to use a fraction of the fitting radius as the initial value
+    for the fitted full-width half-maximum in cases where the moment analysis
+    cannot compute an initial value.
+    (4/15/94 LED)
+
+pkg/images/imarith/iclog.x
+    Changed the mean, median, mode, and zero formats from 6g to 7.5g to
+    insure 5 significant digits regardless of signs and decimal points.
+    (4/13/94, Valdes)
+
+pkg/images/doc/imcombine.hlp
+    Tried again to clarify the scaling as multiplicative and the offseting
+    as additive for file input and for log output.  (3/22/94, Valdes)
+
+pkg/images/imarith/iacclip.gx
+pkg/images/imarith/iccclip.gx
+pkg/images/imarith/iscclip.gx
+    The image sigma was incorrectly computed when an offset scaling is used.
+    (3/8/94, Valdes)
+
+pkg/images/doc/imcombine.hlp
+    The MINMAX example confused low and high.  (3/7/94, Valdes)
+
+pkg/images/geometry/t_geomap.x
+pkg/images/geometry/geofit.x
+pkg/images/geometry/geograph.x
+    Fixed a bug in the geomap code which caused the linear portion of the transformation
+    to be computed incorrectly if the x and y fits had a different functional form.
+    (12/29/93, Davis)
+
+pkg/images/imarith/t_imcombine.x
+pkg/images/imcombine.par
+pkg/images/do/imcombine.hlp
+    The output pixel datatypes now include unsigned short integer.
+    (12/4/93, Valdes)
+
+pkg/images/doc/imcombine.hlp
+    Fixed an error in the example of offseting.  (11/23/93, Valdes)
+
+pkg/images/imfit/fit1d.x
+    When doing operations in place the input and output buffers are the
+    same and the difference and ratio operations assumed they were not
+    causing the final results to be wrong.  (11/16/93, Valdes)
+
+pkg/images/imarith/t_imarith.x
+pkg/images/doc/imarith.hlp
+    If no calculation type is specified then it will be at least real
+    for a division.  Since the output pixel type defaults to the
+    calculation type if not specified this will also result in a
+    real output if dividing two integer images.  (11/12/93, Valdes)
+
+pkg/images/imarith/icgrow.gx
+pkg/images/imarith/icpclip.gx
+pkg/images/imarith/icsclip.gx
+pkg/images/imarith/icaclip.gx
+pkg/images/imarith/iccclip.gx
+pkg/images/imarith/t_imcombine.x
+pkg/images/doc/imcombine.hlp
+    If there were fewer initial pixels than specified by nkeep then the
+    task would attempt to add garbage data to achieve nkeep pixels.  This
+    could occur when using offsets, bad pixel masks, or thresholds.  The
+    code was changed to check against the initial number of pixels rather
+    than the number of images.  Also a negative nkeep is no longer
+    converted to a positive value based on the number of images.  Instead
+    it specifies the maximum number of pixels to reject from the initial
+    set of pixels.  (11/8/93, Valdes)
+
+=======
+V2.10.2
+=======
+
+pkg/images/imarith/icsetout.x
+    Added MWCS calls to update the axis mapping when using the project
+    option in IMCOMBINE.  (10/8/93, Valdes)
+
+pkg$images/imarith/icscale.x
+pkg$images/doc/imcombine.hlp
+    The help indicated that user input scale or zero level factors
+    by an @file or keyword are multiplicative and additive while the
+    task was using then as divisive and subtractive.  This was
+    corrected to agree with the intend of the documentation.
+    Also the factors are no longer normalized.  (9/24/93, Valdes)
+
+pkg$images/imarith/icsetout.x
+    The case in which absolute offsets are specified but the offsets are
+    all the same did not work correctly.  (9/24/93, Valdes)
+
+pkg$images/imfit/imsurfit.h
+pkg$images/imfit/t_imsurfit.x
+pkg$images/imfit/imsurfit.x
+pkg$images/lib/ranges.x
+    Fixed two bugs in the imsurfit task bad pixel rejection code. For low
+    k-sigma rejections factors the bad pixel list could overflow resulting
+    in a segmentation violation or a hung task. Overlapping ranges were
+    not being decoded into a bad pixel list properly resulting in 
+    oscillating bad pixel rejection behavior where certain groups of
+    bad pixels were alternately being included and excluded from the fit.
+    Both bugs are fixed in iraf 2.10.3
+    (9/21/93, Davis)
+
+pkg$images/doc/imcombine.hlp
+    Clarified how bad pixel masks work with the "project" option.
+    (9/13/93, Valdes)
+
+pkg$images/imfit/fit1d.x
+    When the input and output images are the same there was an typo error
+    such that the output was opened separately but then never unmapped
+    resulting in the end of the image not being updated.  (8/6/93, Valdes)
+
+pkg$images/imarith/t_imcombine.x
+    The algorithm for making sure there are enough file descriptors failed
+    to account for the need to reopen the output image header for an
+    update.  Thus when the number of input images + output images + logfile
+    was exactly 60 the task would fail.  The update occurs when the output
+    image is unmapped so the solution was to close the input images first
+    except for the first image whose pointer is used in the new copy of the
+    output image.  (8/4/93, Valdes)
+
+pkg$images/filters/t_mode.x
+pkg$images/filters/t_median.x
+    Fixed a bug in the error trapping code in the median and mode tasks.
+    The call to eprintf contained an extra invalid error code agument.
+    (7/28/93, Davis)
+
+pkg$images/geometry/geomap.par
+pkg$images/geometry/t_geomap.x
+pkg$images/geometry/geogmap.x
+pkg$images/geometry/geofit.x
+    Fixed a bug in the error handling code in geomap which was producing
+    a segmentation violation on exit if the user's coordinate list
+    had fewer than 3 data points. Also improved the error messages
+    presented to the user in both interactive and non-interactive mode.
+    (7/7/93, Davis)
+
+pkg$images/imarith/icgdata.gx
+    There was an indexing error in setting up the ID array when using
+    the grow option.  This caused the CRREJECT/CCDCLIP algorithm to
+    fail with a floating divide by zero error when there were non-zero
+    shifts.  (5/26/93, Valdes)
+
+pkg$images/imarith/icmedian.gx
+    The median calculation is now done so that the original input data
+    is not lost.  This slightly greater inefficiency is required so
+    that an output sigma image may be computed if desired.  (5/10/93, Valdes)
+
+pkg$images/geometry/t_imshift.x
+    Added support for type ushort to the imshift task in cases where the
+    pixel shifts are integral.
+    (5/8/93, Davis)
+
+pkg$images/doc/rotate.hlp
+    Fixed a bug in the rotate task help page which implied that automatic
+    image size computation would occur if ncols or nlines were set no 0
+    instead of ncols and nlines.
+    (4/17/93, Davis)
+
+pkg$images/imarith/imcombine.gx
+    There was no error checking when writing to the output image.  If
+    an error occurred (the example being when an imaccessible imdir was
+    set) obscure messages would result.  Errchks were added.
+    (4/16/93, Valdes)
+
+pkg$images/doc/gauss.hlp
+    Fixed 2 sign errors in the equations in the documentation describing
+    the elliptical gaussian fucntion.
+    (4/13/92, Davis)
+
+pkg/images/imutil/t_imslice.x
+    Removed an error check in the imslice task, which was preventing it from
+    being used to reduce the dimensionality of images where the length of 
+    the slice dimension is 1.0.
+    (2/16/83, Davis)
+
+pkg/images/filters/fmedian.x
+    The fmedian task was printing debugging information under iraf 2.10.2.
+    (1/25/93, Davis)
+
+pkg/images/imarith/icaclip.gx
+pkg/images/imarith/iccclip.gx
+pkg/images/imarith/icpclip.gx
+pkg/images/imarith/icsclip.gx
+    When using mclip=yes and when more pixels are rejected than allowed by
+    the nkeep parameter there was a subtle bug in how the pixels are added
+    back which can result in a segmentation violation.
+	if (nh == n2)  ==>  if (nh == n[i])
+    (1/20/93, Valdes)
+
+
+=======
+V2.10.1
+=======
+
+pkg/images/imarith/t_imcombine.x
+pkg/images/imarith/icaclip.gx
+pkg/images/imarith/iccclip.gx
+pkg/images/imarith/icgrow.gx
+pkg/images/imarith/iclog.x
+pkg/images/imarith/icombine.com
+pkg/images/imarith/icombine.gx
+pkg/images/imarith/icombine.h
+pkg/images/imarith/icpclip.gx
+pkg/images/imarith/icscale.x
+pkg/images/imarith/icsclip.gx
+pkg/images/imarith/icsetout.x
+pkg/images/imcombine.par
+pkg/images/doc/combine.hlp
+    The weighting was changed from using the square root of the exposure time
+    or image statistics to using the values directly.  This corresponds
+    to variance weighting.  Other options for specifying the scaling and
+    weighting factors were added; namely from a file or from a different
+    image header keyword.  The \fInkeep\fR parameter was added to allow
+    controlling the maximum number of pixels to be rejected by the clipping
+    algorithms.  The \fIsnoise\fR parameter was added to include a sensitivity
+    or scale noise component to the noise model.  Errors will now delete
+    the output image.
+    (9/30/92, Valdes)
+
+pkg/images/imutil/imcopy.x
+    Added a call to flush after the status line printout so that the output
+    will appear immediately. (8/19/92, Davis)
+
+pkg/images/filters/mkpkg
+pkg/images/filters/t_fmedian.x
+pkg/images/filters/fmedian.x
+pkg/images/filters/fmd_buf.x
+pkg/images/filters/fmd_maxmin.x
+    The fmedian task could crash with a segmentation violation if mapping
+    was turned off (hmin = zmin and hmax = zmax) and the input image
+    contained data outside the range defined by zmin and zmax. (8/18/92, Davis)
+
+pkg/images/imarith/icaclip.gx
+pkg/images/imarith/iccclip.gx
+pkg/images/imarith/icpclip.gx
+pkg/images/imarith/icsclip.gx
+    There was a very unlikely possibility that if all the input pixels had
+    exactly the same number of rejected pixels the weighted average would
+    be done incorrectly because the dflag would not be set.  (8/11/92, Valdes)
+
+pkg/images/imarith/icmm.gx
+    This procedure failed to set the dflag resulting in the weighted average
+    being computed in correctly.  (8/11/92, Valdes)
+
+pkg/images/imfit/fit1d.x
+    At some point changes were made but not documented dealing with image
+    sections on the input/output.  The changes seem to have left off the
+    final step of opening the output image using the appropriate image
+    sections.  Because of this it is an error to use an image section
+    on an input image when the output image is different; i.e.
+    
+	cl> fit1d dev$pix[200:400,*] junk
+	
+    This has now been fixed.  (8/10/92, Valdes)
+
+pkg/images/imarith/icscales.x
+    The zero levels were incorrectly scaled twice.  (8/10/92, Valdes)
+
+pkg/images/imarith/icstat.gx
+    Contained the statement
+	nv = max (1., (Memi[v2+i] - Memi[v1+i]) / Memi[dv+i] + 1)
+    which is max(real,int).  Changed the 1. to a 1.  (8/10/92, Valdes)
+
+pkg$images/imarith/icaclip.gx
+pkg$images/imarith/iccclip.gx
+pkg$images/imarith/icsclip.gx
+    These files contained multiple cases (ten or so) of constructs such as
+    "max (1., ...)" or "max (0., ...)" where the ... could be either real
+    or double.   In the double cases the DEC compiler complained about a
+    type mismatch since 1. is real.  (8/10/92, Valdes)
+
+pkg$images/imfit/t_imsurfit.x
+    Fixed a bug in the section reading code. Imsurfit is supposed to switch
+    the order of the section delimiters in x and y if x2 < x1 or y2 < 1.
+    Unfortunately the y test was actually "if (y2 < x1)" instead of 
+    "if (y2 < y1)". Whether or not the code actually works correctly
+    depends on the value of x1 relative to y2. This bug was not present
+    in 2.9.1 but is present in subsequent releases. (7/30/92 LED)
+
+=======
+V2.10.1
+=======
+
+pkg$images/filters/t_gauss.x
+    The case theta=90 and ratio > 0.0 but  < 1.0 was producing an incorrect
+    convolution if bilinear=yes, because the major axis sigmas being 
+    input along the x and y axes were sigma and ratio * sigma respectively
+    instead of ratio * sigma and sigma in this case.
+
+pkg$images/imutil/imcopy.x
+    Modified imcopy to write its verbose output to STDOUT instead of
+    STDERR. (6/24/92, Davis)
+
+pkg$images/imarith/imcombine.gx
+    The step where impl1$t is called to check if there is enough memory
+    did not set the return buffer because the values are irrelevant for
+    this check.  However, depending on history, this buffer could have
+    arbitrary values and later when IMIO attempts to flush this buffer,
+    at least in the case of image type coersion, cause arithmetic errors.
+    The fix was to clear the returned buffers.  (4/27/92, Valdes)
+
+pkg$images/imutil/t_imstack.x
+    Modified the imslice task to read the old and write a new axis map.
+    (4/23/92, Davis)
+
+pkg$images/geometry/t_imslice.x
+    Modified the imslice task to read the old and write a new axis map.
+    (4/23/92, Davis)
+
+pkg$images/geometry/t_blkavg.x
+pkg$images/geometry/t_blkrep.x
+    Modified the calls to mw_shift and mw_scale to explicitly set the
+    number of logical axes instead of using the default of 0.
+    (4/23/92, Davis)
+
+pkg$images/geometry/t_imtrans.x
+    Modified imtranspose so that it correctly picks up the axis map
+    and writes it to the output image wcs. (4/23/92, Davis)
+
+pkg$images/register.par
+pkg$images/geotran.par
+pkg$images/doc/register.hlp
+pkg$images/doc/geotran.hlp
+    Changed the default values of the parameters xscale and yscale in
+    the register and geotran tasks from INDEF to 1.0 (4/23/92, Davis)
+
+pkg$images/geometry/t_imtrans.x
+pkg$images/doc/imtranspose.hlp
+    Modified the imtranspose task so it does a true transpose of the
+    axes instead of simply modifying the lterm. (4/8/92, Davis)
+
+pkg$images/iminfo/listpixels.x
+    Added the formats parameter for formatting the output pixel coordinates
+    to the listpixels task. These formats take precedence over the formats
+    stored in the WCS in the image header and the previous default format.
+    (4/7/92, Davis)
+
+pkg$images/imutil/t_imstack.x
+    Added wcs support to the imstack task. (4/2/92, Davis)
+
+pkg$images/iminfo/listpixels.x
+    Modified listpixels so that it will work correctly  if the dimension
+    of the wcs is less than the dimension of the image. (3/16/92, Davis)
+
+pkg$images/geometry/t_geotran.x
+    Modified the rotate, imlintran, register and geotran tasks wcs updating
+    code to deal correclty with dimensionally reduced data. (3/16/92, Davis)
+
+pkg$images/imarith/icalip.gx
+pkg$images/imarith/icclip.gx
+pkg$images/imarith/ipslip.gx
+pkg$images/imarith/icslip.gx
+pkg$images/imarith/icmedian.gx
+    The median calculation with an even number of points for short data
+    could overflow (addition of two short values) and be incorrect.
+    (3/16/92, Valdes)
+
+pkg$images/geometry/t_blkavg.x
+pkg$images/geometry/t_blkrep.x
+    1. Improved the precision of the blkavg task wcs updating code.
+    2. Changed the blkrep task wcs updating code so that it is consistent
+    with blkavg. This means that a blkrep command followed by a blkavg
+    command or vice versa will return the original coordinate system
+    to within machine precision. (3/16/92, Davis)
+
+pkg$images/iminfo/listpixels.x
+    Modified listpixels to print out an error if it could not open the
+    wcs in the image. (3/15/92, Davis)
+
+pkg$images/geometry/t_magnify.x
+    Fixed a bug in the magnify task wcs updating code which was not
+    working correctly for dimensionally reduced images. (3/15/92, Davis)
+
+pkg$images/geometry/t_imtrans.x
+    Fixed a bug in the imtranspose task wcs updating code which was not
+    working correctly for dimensionally reduced images. (3/14/92, Davis)
+
+pkg$images/imarith/icalip.gx
+pkg$images/imarith/icclip.gx
+pkg$images/imarith/icslip.gx
+    There was a bug allowing the number of valid pixels counter to become
+    negative.  Also there was a step which should not be done if the
+    number of valid pixels is less than 1; i.e. all pixels rejected.
+    A test was put in to skip this step.  (3/13/92, Valdes)
+
+pkg$images/iminfo/t_imslice.x
+pkg$images/doc/imslice.hlp
+    Added wcs support to the imslice task.
+    (3/12/92, Davis)
+
+pkg$images/iminfo/t_imstat.x
+    Fixed a bug in the code for computing the standard deviation, kurtosis,
+    and skew, wherein precision was being lost because two of the intermediate
+    variables in the computation were real instead of double precision.
+    (3/10/92, Davis)
+
+pkg$images/iminfo/listpixels.x
+    1. Modified listpixels task to use the MWCS axis "format" attributes 
+    if they are present in the image header.
+    2. Added support for dimensionally reduced images, i.e.
+    images which are sections of larger images and whose coordinate
+    transformations depend on the reduced axes, to the listpixels task. 
+    (3/6/92, Davis)
+
+pkg$images/imarith/t_imcombine.x
+pkg$images/imarith/icsetout.x
+    Changed error messages to say IMCOMBINE instead of ICOMBINE.
+    (3/2/92, Valdes)
+
+pkg$images/imarith/iclog.x
+    Added listing of read noise and gain.  (2/10/92, Valdes)
+
+pkg$images/imarith/icscale.x
+pkg$images/imarith/icpclip.gx
+    1.  Datatype declaration for asumi was incorrect.
+    2.  Reduced the minimum number of images allowed for PCLIP to 3.
+    (1/7/92, Valdes)
+
+pkg$images/imarith/icgrow.gx
+    The first pixel to be checked was incorrectly set to 0 instead of 1
+    resulting in a segvio when using the grow option.  (12/6/91, Valdes)
+
+pkg$images/imarith/icgdata.gx
+pkg$images/imarith/icscale.x
+    Fixed datatype declaration errors found by SPPLINT. (11/22/91, Valdes)
+
+pkg$images/iminfo/t_imstat.x
+    Fixed a bug in the kurtosis computation found by ST.
+    (Davis 10/11/91)
+
+pkg$images/iminfo/t_imstat.x
+pkg$images/doc/imstat.hlp
+    Corrected a bug in the mode computation in imstatistics. The parabolic
+    interpolation correction for computing the histogram peak was being
+    applied in the wrong direction. Note that for dev$pix the wrong answer
+    is actually closer to the expected answer than the correct answer
+    due to binning effects. 
+    (Davis 9/24/91)
+
+pkg$images/filters/t_gauss.x
+    The code which computes the gaussian kernel was producing a divide by
+    zero error if ratio=0.0 and bilinear=yes (2.10 version only).
+    (Davis 9/18/91)
+
+pkg$images/doc/magnify.hlp
+    Corrected a bug in the magnify help page.
+    (Davis 9/18/91)
+
+pkg$images/imarith/icsclip.gx
+pkg$images/imarith/icaclip.gx
+pkg$images/imarith/iccclip.gx
+    There was a typo,  Memr[d[k]+k] --> Memr[d[j]+k].  (9/17/91, Valdes)
+
+pkg$images/imarith/icstat.gx
+pkg$images/imarith/icmask.x
+    The offsets were used improperly in computing image statistics.
+    (Valdes, 9/17/91)
+
+pkg$images/geometry/t_imshift.x
+    The shifts file pointer was not being correctly initialized to NULL 
+    in the case where no shifts file was declared. When the task
+    was invoked repeatedly from a script, this could result in an array being
+    referenced, for which space had not been previously allocated.
+    (Davis 7/29/91)
+
+pkg$images/imarith/imc* -
+pkg$images/imarith/ic* +
+pkg$images/imarith/t_imcombine.x
+pkg$images/imarith/mkpkg
+pkg$images/imarith/generic/mkpkg
+pkg$images/imcombine.par
+pkg$images/doc/imcombine.hlp
+    Replaced old version of IMCOMBINE with new version supporting masks,
+    offsets, and new algorithms.  (Valdes 7/19/91)
+
+pkg$images/iminfo/imhistogram.x
+    Imhistogram has been modified to print the value of the middle of
+    histogram bin instead of the left edge if the output type is list 
+    instead of plot. (Davis 6/11/91)
+
+pkg$images/t_imsurfit.x
+    Modified the sections file reading code to check the order of the
+    x1 x2 y1 y2 parameters and switch (x1,x2) or (y1,y2) if x2 < x1 or
+    y2 < y1 respectively. (Davis 5/28/91)
+
+pkg$images/listpixels.par
+pkg$images/iminfo/listpixels.x
+pkg$images/doc/listpixels.hlp
+    Modified the listpixels task to be able to print the pixel coordinates
+    in logical, physical or world coordinates. The default coordinate
+    system is still logical as before. (Davis 5/17/91)
+
+pkg$images/images.par
+pkg$images/doc/minmax.hlp
+pkg$images/imutil/t_minmax.x
+pkg$images/imutil/minmax.x
+    Minmax was modified to do the minimum and maximum values computations
+    in double precision or complex instead of real if the input image
+    pixel type is double precision or complex. Note that the minimum and
+    maximum header values are still stored as real however.
+    (Davis 5/16/91)
+
+imarith/t_imarith.x
+    There was a missing statement to set the error flag if the image
+    dimensions did not  match.  (5/14/91, Valdes)
+
+doc/imarith.hlp
+    Fixed some formatting problems in the imarith help page. (5/2/91 Davis)
+
+imarith$imcombine.x
+    Changed the order in which images are unmapped to have the output images
+    closed last.  This is to allow file descriptors for the temporary image
+    used when updating STF headers.  (4/22/91, Valdes)
+
+pkg$images/geometry/t_blkavg.x
+pkg$images/geometry/blkavg.gx
+pkg$images/geometry/blkavg.x
+    The blkavg task was partially modified to support complex image data. 
+    The full modifications cannot be made because of an error in abavx.x
+    and the missing routine absux.x.
+    (4/18/91 Davis)
+
+pkg$images/geometry/geofit.x
+    The x and y fits cross-terms switch was not being set correctly to "yes"
+    in the case where xxorder=2 and xyorder=2 or in the case where yxorder=2
+    and yyorder=2.
+    (4/9/91 Davis)
+
+pkg$images/geometry/geogmap.x
+    Modified the line which prints the geometric parameters to use the
+    variable name xshift and yshift instead of delx and dely.
+    (4/9/91 Davis)
+
+pkg$images/imfit/imsurfit.x
+    Fixed a bug in the pixel rejection code which occurred when upper was >
+    0.0 and lower = 0.0 or lower > 0 and upper = 0.0. The problem was that
+    the code was simply setting the rejection limits to the computed sigma
+    times the upper and lower parameters without checking for the 0.0
+    condition first.  In the first case this results in all points with
+    negative residuals being rejected and in the latter all points with
+    positive residuals are rejected.
+    (2/25/91 Davis)
+
+pkg$images/doc/hedit.hlp
+pkg$images/doc/hselect.hlp
+pkg$images/doc/imheader.hlp
+pkg$images/doc/imgets.hlp
+    Added a reference to imgets in the SEE ALSO sections of the hedit and
+    hselect tasks.
+    Added a reference to hselect and hedit in the SEE ALSO sections of the
+    imheader and imgets tasks.
+    (2/22/91 Davis)
+
+pkg$images/gradient.hlp
+pkg$images/laplace.hlp
+pkg$images/gauss.hlp
+pkg$images/convolve.hlp
+pkg$images/gradient.par
+pkg$images/laplace.par
+pkg$images/gauss.par
+pkg$images/convolve.par
+pkg$images/t_gradient.x
+pkg$images/t_laplace.x
+pkg$images/t_gauss.x
+pkg$images/t_convolve.x
+pkg$images/convolve.x
+pkg$images/xyconvolve.x
+pkg$images/radcnv.x
+    The convolution operators were modified to run more efficiently in
+    certain cases. The LAPLACE task was modified to make use of the
+    radial symmetry of the convolution kernel in the y direction as well
+    as the x direction resulting in a modest speedup in execution time.
+    A new parameter bilinear was added to the GAUSS and CONVOLVE tasks.
+    By default and if appropriate mathematically,  GAUSS now makes use of
+    the bilinearity or separability of the Gaussian function,
+    to separate the 2D convolution in x and y into two equivalent
+    1D convolutions in x and y, resulting in a considerable speedup
+    in execution time. Similarly the user can know program CONVOLVE to
+    compute a bilinear convolution instead of a full 2D 1 if appropriate.
+    (1/29/91 Davis)
+
+pkg$images/filters/t_convolve.x
+    CONVOLVE was not decoding the legal 1D kernel "1.0 2.0 1.0" correctly
+    although the alternate form "1.0 2.0 1.0;" worked. Leading
+    blanks in string kernels as in for example " 1.0 2.0 1.0" also generated
+    and error.  Fixed these bugs and added some  additional error checking code.
+    (11/28/90 Davis)
+
+pkg$images/doc/gauss.hlp
+    Added a detailed mathematical description of the gaussian kernel used
+    by the GAUSS task to the help page.
+
+pkg$images/images.hd
+pkg$images/rotate.cl
+pkg$images/imlintran.cl
+pkg$images/register.cl
+pkg$images/register.par
+	Added src="script file name" entries to the IMAGES help database
+	for the tasks ROTATE, IMLINTRAN, and REGISTER. Changed the CL
+	script for REGISTER to a procedure script to remove the ugly
+	local variable declarations. Added a few comments to the scripts.
+	(12/11/90, Davis)
+
+pkg$images/iminfo/imhistogram.x
+    Added a new parameter binwidth to imhistogram. If binwidth is defined
+    it determines the histogram resolution in intensity units, otherwise
+    nbins determines the resolution as before. (10/26/90, Davis)
+
+pkg$images/doc/sections.hlp
+    Clarified what is meant by an image template and that the task itself
+    does not check whether the specified names are actually images.
+    The examples were improved.  (10/3/90, Valdes)
+
+pkg$images/doc/fit1d.hlp
+    Changed lines to columns in example 2.  (10/3/90, Valdes)
+
+pkg$images/imarith/imcscales.x
+    When an error occured while parsing the mode section the untrapped error
+    caused further problems downstream.  Because it would require adding
+    lots of errchks to cause the program to gracefully abort I instead made
+    it a warning.  (10/2/90, Valdes)
+
+pkg$images/imutil/hedit.x
+    Hedit was computing but not using min_lenarea. If the user specified
+    a min_lenuserarea greater than the default of 28800 then the default
+    was being used instead of the larger number.
+
+pkg$imarith/imasub.gx
+    The case of subtracting an image from the constant zero had a bug
+    which is now fixed.  (8/14/90, Valdes)
+
+pkg$images/t_imtrans.x
+    Modified the imtranspose task so it will work on type ushort images.
+    (6/6/90 Davis)
+
+pkg$images
+    Added world coordinate system support to the following tasks: imshift,
+    shiftlines, magnify, imtranspose, blkrep, blkavg, rotate, imlintran,
+    register and geotran. The only limitation is that register and geotran
+    will only support simple linear transformations. 
+    (2/24/90 Davis)
+    
+pkg$images/geometry/geotimtran.x
+    Fixed a problem in the boundary extension "reflect" option code for small
+    images which was causing odd values to be inserted at the edges of the
+    image.
+    (2/14/90 Davis)
+
+pkg$images/iminfo/imhistogram.x
+    A new parameter "hist_type" was added to the imhistogram task giving
+    the user the option of plotting the integral, first derivative and
+    second derivative of the histogram as well as the normal histogram.
+    Code was contributed by Rob Seaman.
+    (2/2/90 Davis)
+
+pkg$images/geometry/geogmap.x
+    The path name of the help file was being erroneously renamed with
+    the result that when users ran the double precision version of the
+    code they could not find the help file.
+    (26/1/90 Davis)
+
+pkg$images/filters/t_boxcar.x,t_convolve.x
+    Added some checks for 1-D images.
+    (1/20/90 Davis)
+
+pkg$images/iminfo/t_imstat.x,imstat.h
+    Made several minor bug fixes and alterations in the imstatistics task
+    in response to user complaints and suggestions.
+
+    1. Changed the verbose parameter to the format parameter. If format is
+    "yes" (the default) then the selected fields are printed in fixed format
+    with column labels. Other wise the fields are printed in free format
+    separated by 2 blanks. This fixes the problem of fields running together.
+
+    2. Fixed a bug in the code which estimates the median from the image
+    histogram by linearly interpolating around the midpt of the integrated
+    histogram. The bug occurred when more than half the pixels were in the
+    first bin.
+
+    3. Added a check to ensure that the number of fields did not overflow
+    the fields array.
+
+    4. Removed the extraneous blank line printed after the title.
+
+    5. The pound sign is now printed at the beginning of the column header
+    string regardless of which field is printed first. In the previous
+    versions it was only being printed if the image name field was
+    printed first.
+
+    6. Changed the name of the median field to midpt in response to user
+    confusions about how the median is computed.
+
+    (1/20/90, Davis)
+
+pkg$images/imutil/t_imslice.hlp
+    The imslice was not correctly computing the number of lines in the
+    output image in the case where the slice dimension was 1.
+    (12/4/89, Davis)
+
+pkg$images/doc/imcombine.hlp
+    Clarified and documented definitions of the scale, offset, and weights.
+    (11/30/89, Valdes)
+
+pkg$images/geometry/geotran.x
+    High order surfaces of a certain functional form could occasionally
+    produce out of bounds pixel errors. The bug was caused by not properly
+    computing the distortion of the image boundary for higher order
+    surfaces.
+    (11/21/89, Davis)
+
+pkg$images/geometry/imshift.x
+    The circulating buffer space was not being freed after each execution
+    of IMSHIFT. This did not cause an error in execution but for a long
+    list of frames could result in alot of memory being tied up.
+    (10/25/89, Davis)
+
+pkg$images/imarith/t_imarith.x
+    IMARITH is not prepared to deal with images sections in the output.
+    It used to look for '[' to decide if the output specification included
+    and image section.  This has been changed to call the IMIO procedure
+    imgsection and check if a non-null section string is returned.
+    Thus it is up to IMIO to decide what part of the image name is
+    an image section.  (9/5/89, Valdes)
+
+pkg$images/imarith/imcmode.gx
+    Fixed bug causing infinite loop when computing mode of constant value
+    section.  (8/14/89, Valdes)
+
+====
+V2.8
+====
+
+pkg$images/iminfo/t_imstat.x
+    Davis, Jun 15, 1989
+    Added a couple of switches to that skew and kurtosis are not computed
+    if they are not to be printed.
+
+pkg$images/iminfo/t_imstat.x
+    Davis, Jun 14, 1989
+    A simple mod was made to the skew and kurtosis computation to avoid
+    divide by zero errors in case of underflow.
+
+pkg$images/imutil/chpixtype.par
+    Davis, Jun 13, 1989
+    The parameter file has been modified to accept an output pixel
+    type of ushort.
+
+pkg$images/imarith/imcombine.gx
+    Valdes, Jun 2, 1989
+    A new scheme to detect file errors is now used.
+
+pkg$images/imfit/t_imsurfit.x
+    Davis, Jun 1, 1989
+    1. If the user set regions = "sections" but the sections file
+    did not exist the task would go into an infinite loop. The problem
+    was a missing error check on the open statement.
+
+pkg$images/iminfo/imhistogram.x,imhistogram.par
+    Davis, May 31, 1989
+    A new version of imhistogram has been installed. These mods have
+    been made over a period of a month by Doug Tody and Rob Seaman.
+    The mods include
+    1. An option to turn off log scaling of the y axis of the histogram plot.
+    2. A new autoscale parameter which avoids aliasing problems for integer
+       data.
+    3. A new parameter top_close which resolves the ambiguity in the top
+       bin of the histogram.
+
+pkg$images/imarith/imcombine.gx
+   Valdes, May 9, 1989
+   Because a file descriptor was not reserved for string buffer operations
+   and a call to stropen in cnvdate was not error checked the task would
+   hang when more than 115 images were combined.  Better error checking
+   was added and now an error message is printed when the maximum number
+   of images that can be combined is exceeded.
+
+pkg$images/imarith/t_imarith.x
+   Valdes, May 6, 1989
+   Operations in which the output image has an image section are now
+   skipped with a warning message.
+
+pkg$images/imarith/sigma.gx
+pkg$images/imarith/imcmode.gx
+    Valdes, May 6, 1989
+    1.  The weighted sigma was being computed incorrectly.
+    2.  The argument declarations were wrong for integer input images.
+	Namely the mean vector is always real.
+    3.  Minor change to imcmode.gx to return correct datatype.
+
+pkg$images/imstack,imslice
+    Davis, April 1, 1989
+    The proto images tasks imstack and imslice have been moved from the
+    proto package to the images package. Imstack is unchanged except that
+    it now supports the image data types USHORT and COMPLEX. Imslice has
+    been modified to allow slicing along any dimension of the image instead
+    of just the highest dimension.
+
+pkg$images/imstatistics.
+    Davis, Mar 31, 1989
+    1. A totally new version of the imstatistics task has been written
+    and replaces the old version. The new task allows the user to select
+    which statistical parameters to compute and print. These include
+    the mean, median, mode, standard deviation, skew, kurtosis and the
+    minimum and maximum pixel values.
+
+pkg$images/imhistogram.par
+pkg$images/iminfo/imhistogram.x
+pkg$images/doc/imhistogram.hlp
+    Davis, Mar 31, 1989
+    1. The imhistogram task has been modified to plot "box" style histograms
+    as well as "line" type histograms. Type "line" remains the default.
+
+pkg$images/geometry/geotran.par,register.par,geomap.par
+pkg$images/doc/geomap.hlp,register.hlp,geotran.hlp
+    Davis, Mar 6, 1989
+    1. Improved the parameter prompting in GEOMAP, REGISTER and GEOTRAN
+    and improved the help pages.
+    2. Changed GEOMAP database quantities "xscale" and "yscale" to "xmag"
+    and "ymag" for consistency . Geotran was changed appropriately.
+
+pkg$images/imarith/imcmode.gx
+    For short data a short variable was wraping around when there were
+    a significant number of saturated pixels leading to an infinite loop.
+    The variables were made real regardless of the image datatype.
+    (3/1/89, Valdes)
+ 
+pkg$images/imutil/imcopy.x
+    Davis, Feb 28, 1989
+    1. Added support for type USHORT to the imcopy task. This is a merged
+    ST modification.
+
+pkg$images/imarith/imcthreshold.gx
+pkg$images/imcombine.par
+pkg$images/doc/imcombine.hlp
+pkg$images/imarith/imcscales.x
+    Valdes, Feb 16, 1989
+    1.  Added provision for blank value when all pixels are rejected by the
+        threshold.
+    2.  Fixed a bug that was improperly scaling images in the threshold option.
+    3.  The offset printed in the log now has the opposite sign so that it
+	is the value "added" to bring images to a common level.
+
+pkg$images/imfit/imsurfit.x
+    Davis, Feb 23, 1989
+    Fixed a bug in the median fitting code which could cause the porgram
+    to go into an infinite loop if the region to be fitted was less than
+    the size of the whole image.
+
+pkg$images/geometry/t_magnify.x
+    Davis, Feb 16, 1989
+    Modified magnify to work on 1D images as well as 2D images. The
+    documentation has been updated.
+
+pkg$images/geometry/t_geotran.x
+    Davis, Feb 15, 1989
+    Modified the GEOTRAN and REGISTER tasks so that they can handle a list
+    of transform records one for each input image.
+
+pkg$images/imarith/imcmode.gx
+    Valdes, Feb 8, 1989
+    Added test for nx=1.
+ 
+pkg$images/imarith/t_imcombine.x
+    Valdes, Feb 3, 1989
+    The test for the datatype of the output sigma image was wrong.
+ 
+pkg$images/iminfo/listpixels.x,listpixels.par
+    Davis, Feb 6, 1989
+    The listpixels task has been modified to print out the pixels for a 
+    list of images instead of a single image only. A title line for each
+    image listed can optionally be printed on the standard output if
+    the new parameter verbose is set to yes.
+
+pkg$images/geometry/t_imshift.x
+    Davis, Feb 2, 1989
+    Added a new parameter shifts_file to the imshift task. Shifts_file
+    is the name of a text file containing the the x and yshifts for
+    each input image to be shifted. The number of input shifts must
+    equal the number of input images.
+
+pkg$images/geometry/t_geomap.x
+    Davis, Jan 17, 1989
+    Added an error message for the case where the coordinates is empty
+    of there are no points in the specified data range. Previously the
+    task would proceed to the next coordinate file without any message.
+
+pkg$images/geometry/t_magnify.x
+    Davis, Jan 14, 1989
+    Added the parameter flux conserve to the magnify task to bring it into
+    line with all the other geometric transformation tasks.
+
+pgk$images/geometry/geotran.x,geotimtran.x
+    Davis, Jan 2, 1989
+    A bug was fixed in the flux conserve code. If the x and y reference
+    coordinates  are not in pixel units and are not 1 then
+    the computed flux per pixel was too small by xscale * yscale.
+
+pkg$images/filters/acnvrr.x,convolve.x,boxcar.x,aboxcar.x
+    Davis, Dec 27, 1988
+    I changed the name of the acnvrr procedure to cnv_radcnvr to avoid
+    a name conflict with a vops library procedure. This only showed
+    up when shared libraries were implemented. I also changed the name
+    of the aboxcarr procedure to cnv_aboxr to avoid conflict with the
+    vops naming conventions.
+
+pkg$images/imarith/imcaverage.gx
+    Davis, Dec 22, 1988
+    Added an errchk statement for imc_scales and imgnl$t to stop the
+    program bombing with segmentation violations when mode <= 0.
+
+pkg$images/imarith/imcscales.x
+    Valdes, Dec 8, 1988
+    1.  IMCOMBINE now prints the scale as a multiplicative quantity.
+    2.  The combined exposure time was not being scaled by the scaling
+	factors resulting in a final exposure time inconsistent with the
+	data.
+
+pkg$images/iminfo/imhistogram.x
+    Davis, Nov 30, 1988
+    Changed the list+ mode so that bin value and count are printed out instead
+    of bin count and value. This makes the plot and list modes compatable.
+
+pkg$images/iminfo/t_imstat.x
+    Davis, Nov 17, 1988
+    Added the n=n+1 back into the inner loop of imstat.
+
+pkg$images/geotran.par,register.par
+    Davis, Nov 11 , 1988
+    Fixed to glaring errors in the parameter files for register and geotran.
+    Xscale and yscale were described as pixels per reference unit when
+    they should be reference units per pixel. The appropriate bug fix has been
+    made.
+
+pkg$images/geometry/t_geotran.x
+    Davis, November 7, 1988
+    The routine gsrestore was not being error checked. If either of the
+    input x or y coordinate surface was linear and the other was not,
+    the message came back GSRESTORE: Illegal x coordinate. This bug has
+    been fixed.
+
+pkg$images/imarith/imcombine.gx
+    Valdes, October 19, 1988
+    A vops clear routine was not called generically causing a crash with
+    double images.
+
+pkg$images/filters/t_fmedian.x,t_median.x,t_fmode.x,t_mode.x,t_gradient.x
+    t_gauss.x,t_boxcar.x,t_convolve.x,t_laplace.x
+    Davis, October 4, 1988
+    I fixed a bug in the error handling code for the filters tasks. If
+    and error occurred during task execution and the input image name was
+    the same as the output image name then the input image was trashed.
+
+pkg$images/imarith/imcscales.gx
+    Valdes, September 28, 1988
+    It is now an error for the mode to be nonpositive when scaling or weighting.
+
+pkg$images/imarith/imcmedian.gx
+    Valdes, August 16, 1988
+    The median option was selecting the n/2 value instead of (n+1)/2.  Thus,
+    for an odd number of images the wrong value was being determined for the
+    median.
+
+pkg$images/geometry/t_imshift.x
+    Davis, August 11, 1988
+    1. Imshift has been modified to uses the optimized code if nearest
+    neighbour interpolation is requested. A nint is done on the shifts
+    before calling the quick shift routine.
+    2. If the requested pixel shift is too large imshift will now
+    clean up any pixelless header files before continuing execution.
+
+pkg$images/geometry/blkavg.gx
+    Davis, July 13, 1988
+    Blkavg has been fixed so that it will work on 1D images.
+
+pkg$images/geometry/t_imtrans.x,imtrans.x
+    Davis, July 12, 1988
+    Imtranspose has been modified to work on complex images.
+
+pkg$images/imutil/t_chpix.x
+    Davis, June 29, 1988
+    A new task chpixtype has been added to the images package. Chpixtype
+    changes the pixel types of a list of images to a specified output pixel
+    type. Seven data types are supported "short", "ushort", "int", "long"
+    "real" and "double".
+
+pkg$images/geometry/rotate.cl,imlintran.cl,t_geotran.x
+    Davis, June 10, 1988
+    The rotate and imlintran scripts have been rewritten to use procedure
+    scripts. This removes all the annoying temporary cl variables which
+    appear when the user does an lpar. In previous versions of these
+    two tasks the output was restricted to being the same size as the input
+    image. This is still the default case, but the user can now set the
+    ncols and nrows parameters to the desired output size. I ncols or nlines
+    < 0 then then the task compute the output image size required to contain
+    the whole input image.
+
+pkg$images/filters/t_convolve.x,t_laplace.x,t_gradient.x,t_gauss.x,convolve.x
+    Davis, June 1, 1988
+    The convolution operators laplace, gauss and convolve have been modified
+    to make use of radial symmetry in the convolution kernel. In gauss and
+    laplace the change is transparent to the user. For the convolve operator
+    the user must indicate that the kernel is radially symmetric by setting
+    the parameter radsym. For kernels of 7 by 7 or greater the speedup
+    in timings is on the order of 30% on the Vax 750 with the fpa.
+
+pkg$images/imarith/imcmode.gx
+    Valdes, Apr 11, 1988
+    1.  The use of a mode sections was handled incorrectly.
+
+pkg$images/imfit/fit1d.x
+    Valdes, Jan 4, 1988
+    1.  Added an error check for a failure in IMMAP.  The missing error check
+	caused FIT1D to hang when a bad input image was specified.
+
+pkg$images/magnify.par
+pkg$images/imcombine.par
+pkg$images/imarith/imcmode.gx
+pkg$images/doc/imarith.hlp
+    Valdes, Dec 7, 1987
+    1.  Added option list to parameter prompts.
+    2.  Fixed minor typo in help page
+    3.  The mode calculation in IMCOMBINE would go into an infinite loop
+	if all the pixel values were the same.  If all the pixels are the
+	same them it skips searching for the mode and returns the constant
+	number.
+
+pkg$images/geometry/geotimtran.x
+    Davis, Nov 25, 1987
+    1. A bug in the boundary extension = wrap option was found in the
+    IMLINTRAN task. The problem occured in computing values for out of
+    bounds pixels in the range 0.0 < x < 1.0, ncols < x < ncols + 1.0,
+    0.0 < y < 1.0 and nlines < y < nlines + 1.  The computed coordinates
+    were falling outside the boundaries of the interpolation array.
+
+pkg$images/geometry/t_geomap.x,geograph.x
+    Davis, Nov 19, 1987
+    1. The geomap task now writes the name of the output file into the database.
+    2. Rotation angles of 360. degrees have been altered to 0 degrees.
+
+pkg$images/imfit/t_imsurfit.x,imsurfit.x
+pkg$images/lib/ranges.x
+    Davis, Nov 2, 1987
+    A bug in the regions fitting option of the IMSURFIT task has been found
+    and fixed. This bug would occur when the user set the regions parameter
+    to sections and then listed section which overlapped each other. The
+    modified ranges package was not handling the overlap correctly and
+    computing a number of points which was incorrect.
+
+pkg$images/imarith/* +
+    Valdes, Sep 30, 1987
+    The directory was reorganized to put generic code in the subdirectory
+    generic.
+
+    A new task called IMCOMBINE has been added.  It provides for combining
+    images by a number of algorithms, statistically weighting the images
+    when averaging, scaling or offsetting the images by the exposure time
+    or image mode before combining, and rejecting deviant pixels.  It is
+    almost fully generic including complex images and works on images of
+    any dimension.
+
+pkg$images/geometry/geotran.x
+    Davis, Sept 3, 1987
+    A bug in the flux conserving algorithm was found in the geotran code.
+    The symptom was that the flux of the output image occasionally was
+    negative. This would happen when two conditions were met, the transformation
+    was of higher order than a simple rotation, magnification, translation
+    and an axis flip was involved. The mathematical interpretation of this
+    bug is that the coordinate surface had turned upside down. The solution
+    for people running systems with this bug is to multiply there images
+    by -1.
+
+pkg$images/imfit/imsurfit.h,t_imsurfit.x
+    Davis, Aug 6, 1987
+    A new option was added to the parameter regions in the imsurfit task.
+    Imsurfit will now fit a surface to a single circular region defined
+    by an x and y center and a radius.
+
+pkg$images/geometry/geotimtran.x
+    Davis, Jun 15, 1987
+    Geotran and register were failing when the output image number of rows
+    and columns was different from the input number of rows and columns.
+    Geotran was mistakenly using the input images sizes to determine the
+    number of output lines that should be produced. The same problem occurred
+    when the values of the boundary pixels were being computed. The program
+    was using the output image dimensions to compute the boundary pixels
+    instead of the input image dimensions.
+
+pkg$images/geometry/geofit.x,geogmap.x
+    Davis, Jun 11, 1987
+    A bug in the error checking code in the geomap task was fixed. The
+    condition of too few points for a reasonable was not being trapped
+    correctly. The appropriate errchk statements were added.
+
+pkg$images/geomap.par
+    Davis, Jun 10, 1987
+    The default fitting function was changed to polynomial. This will satisfy
+    most users who wish to do shifts, rotations, and magnifications and
+    avoid the neccessity of correctly setting the xmin, xmax, ymin, and ymax
+    parameters. For the chebyshev and legendre polynomial functions these
+    parameters must be explicitly set.  For reference coordinates in pixel
+    units the normal settings are 1, ncols, 1 and nlines respectively.
+
+pkg$images/iminfo/hselect.x,imheader.x,images$/imutil/hselect.x
+    Davis, Jun 8, 1987
+    Imheader has been modified to open an image with the default min_lenuserarea
+    Hselect and hedit will now open the image setting the user area to the
+    maximum of 28800 chars or the min_lenuser environment variable.
+
+pkg$images/iminfo/t_imstat.x
+    Davis, May 22, 1987
+    An error in the image minimum computation was corrected. This error
+    would show up most noiticeably if imstat was run on a 1 pixel image.
+    The min value would be left set to MAX_REAL.
+
+pkg$images/filters/mkpkg
+    Davis, May 22, 1987
+    I added mach.h to the dependency file list of t_fmedian.x and
+    recompiled. The segmentation violations I had been getting in the
+    program disappeared.
+
+pkg$images/t_shiftlines.x,shiftlines.x
+    Davis, April 15, 1987
+    1. I changed the names of the procedures shiftlines and shiftlinesi
+    to sh_lines and sh_linesi. When the original names were contracted
+    to 6 letter fortran names they became shifti and shifts which just
+    so happens to collide with shifti and shifts in the subdirectory
+    osb. On VMS this was causing problems with the shareable libraries.
+    If images was linked with -z there was no problem.
+
+pkg$images/imarith/t_imsum.x
+    Valdes, March 24, 1987
+    1.  IMSUM was failing to unmap images opened to check image dimensions
+	in a quick first pass through the image list.  This is probably
+	the source of the out of files problem with STF images.  It may
+	be the source of the out of memory problem reported from AOS/IRAF.
+
+pkg$images/imfit/fit1d.x
+pkg$images/imfit/mkpkg
+    Valdes, March 17, 1987
+    1. Added error checking for the illegal operation in which both input
+       and output image had an image section.  This was causing the task
+       to crash.  The task now behaves properly in this circumstance and
+       even allows the fitted output to be placed in an image section of
+       an existing output image (even different than the input image
+       section) provided the input and output images have the same sizes.
+
+pkg$images/t_convolve.x
+    Davis, March 3, 1987
+    1. Fixed the kernel decoding routine in the convolve task so that
+    it now recognizes the row delimter character in string entry mode.
+
+pkg$images/geometry,filters
+    Davis, February 27, 1987
+    1. Changed all the imseti (im, TY_BNDRYPIXVAL, value) calls to imsetr.
+
+pkg$images/t_minmax.x,minmax.x
+    Davis, February 24, 1987
+    1. Minmax has been changed to compute the minimum and maximum pixel
+    as well as the minimum and maximum pixel values. The pixels are output
+    in section notation and stored in the minmax parameter file.
+
+pkg$images/t_magnify.x
+    Davis, February 19, 1987
+    1. Magnify was aborting with the error MSIFIT: Too few datapoints
+    when trying to reduce an image using the higher order interpolants
+    poly3, poly5 and spline3. I increased the NEDGE defined constant
+    from 2 to three and modified the code to use the out of bounds
+    imio.
+
+pkg$images/geograph.x,geogmap.x
+    Davis, February 17, 1987
+    1. Geomap now uses the gpagefile routine to page the .keys file.
+    The :show command deactivates the workstation before printing a
+    block of text and reactivates it when it is finished.
+
+pkg$images/geometry/geomap,geotran
+    Davis, January 26, 1987
+    1. There have been substantial changes to the geomap, and geotrans
+    tasks and those tasks rotate, imlintran and register which depend
+    on them.
+    2. Geomap has been changed to be able to compute a transformation
+    in both single and double precision.
+    3. The geotran code has been speeded up considerably. A simple rotate
+    now takes 70 seconds instead of 155 seconds using bilinear interpolation.
+    4. Two new cl parameters nxblock and nyblock have been added to the
+    rotate, imlintran, register and geotran tasks. If the output image
+    is smaller than these parameters then the entire output image
+    is computed at once. Otherwise the output image is computed in blocks
+    nxblock by nyblock in size.
+    5. The 3 geotran parameters rotation, scangle and flip have been replaced
+    with two parameters xrotation and yrotation which serve the same purpose.
+
+pkg$images/geometry/t_shiftlines.x,shiftlines.x
+    Davis, January 19, 1987
+    1. The shiftlines task has been completely rewritten. The following
+    are the major changes.
+    2. Shiftlines now makes use of the imio boundary extension operations.
+    Therefore the four options: nearest pixel, reflect, wrap and constant
+    boundary extension are available.
+    3. The interpolation code has been vectorised. The previous version
+    was using the function call asieval for every output pixel evaluated.
+    The asieval call were replaced with asivector calls.
+    4. An extra CL parameter constant to support constant boundary
+    exension was added.
+    5. The shiftlines help page was modified and the date changed to
+    January 1987.
+
+pkg$images/imfit/imsurfit.x
+    Davis, January 12, 1987
+    1. I changed the amedr call to asokr calls. For my application it did
+    not matter whether the input array is left partially sorted and the asokr
+    routine is more efficient.
+
+pkg$images/lib/pixlist.x
+    Davis, December 12, 1986
+    1. A bug in the pl_get_ranges routine caused the routine to fail when the
+    number of ranges got too large. The program could not detect the end of
+    the ranges and would go into an infinite loop.
+
+pkg$images/iminfo/t_imstat.x
+    Davis, December 3, 1986
+    1. Imstat was failing on constant images because finite machine precision
+    could result in a negative sigma squared. Added a check for this condition.
+
+pkg$images/filters/fmode.x
+    Davis, October 27, 1986
+    1. Added a check for 0 data range before calling amapr.
+
+pkg$images/imarith/imsum.gx
+    Valdes, October 20, 1986
+    1. Found and fixed bug in this routine which caused pixel rejection
+    to fail some fraction of the time.
+
+pkg$images/geometry/blkrp.gx
+    Valdes, October 13, 1986
+    1.  There was a bug when the replication factor for axis 1 was 1.
+
+pkg$images/iminfo/imhistogram.x
+    Hammond, October 8, 1986
+    1. Running imhistogram on a constant valued image would result in
+       a "floating divide by zero fault" in ahgm.  This condition is
+       now trapped and a warning printed if there is no range in the data.
+
+pkg$images/tv/doc/cvl.hlp
+    Valdes, October 7, 1986
+    1.  Typo in V2.3 documentation fixed: "zcale" -> "zscale".
+
+pkg$images/fit1d.par
+    Valdes, October 7, 1986
+    1.  When querying for the output type the query was:
+
+Type of output (fit, difference, ratio) (fit|difference|ratio) ():
+
+    	The enumerated values were removed since they are given in the
+	prompt string.
+
+pkg$images/imarith/t_imsum.x
+pkg$images/imarith/imsum.gx
+pkg$images/do/imsum.hlp
+    Valdes, October 7, 1986
+    1.  Medians or pixel rejection with more than 15 images is now
+	correct.  There was an error in buffering.
+    2.  Averages of integer datatype images are now correct.  The error
+	was caused by summing the pixel values divided by the number
+	of images instead of summing the pixel values and then dividing
+	by the number of images.
+    3.  Option keywords may now be abbreviated.
+    4.  The output pixel datatype now defaults to the calculation datatype
+	as is done in IMARITH.  The help page was modified to indicate this.
+    5.  Dynamic memory is now used throughout to reduce the size of the
+	executable.
+    6.  The bugs 1-2 are present in V2.3 and not in V2.2.
+
+pkg$images/imarith/t_imarith.x
+pkg$images/imarith.par
+pkg$images/doc/imarith.hlp
+    Valdes, October 6, 1986
+    1.  The parameter "debug" was changed to "noact".  "debug" is reserved
+	for debugging information.
+    2.  The output pixel type now defaults to the calculation datatype.
+    3.  The datatype of constant operands is determined with LEXNUM.  This
+	fixes a bug in which a constant such as "1." was classified as an
+	integer.
+    4.  Trailing whitespace in the string for a constant operand is allowed.
+	This fixes a bug with using "@" files created with the task FIELDS
+	from a table of numbers.  Trailing whitespace in image names is
+	not checked for since this should be taken care of by lower level
+	system services.
+    5.  The reported bug with the "max" operation not creating a pixel file
+	was the result of the previous round of changes.  This has been
+	corrected.  This problem does not exist in the released version.
+    6.  All strings are now dynamically allocated.  Also IMTOPENP is used
+	to open a CL list directly.
+    7.  The help page was revised for points (1) and (2).
+
+pkg$images/fmode.par
+pkg$images/fmd_buf.x
+pkg$images/med_sort.x
+    Davis, September 29, 1986
+    1. Changed the default value of the unmap parameter in fmode to yes. The
+       documentation was changed and the date modified.
+    2. Added a test to make sure that the input image was not a constant
+       image in fmode and fmedian.
+    3. Fixed the recently added swap macro in the sort routines which
+       was giving erroneous results for small boxes in tasks median and mode.
+
+pkg$images/imfit/fit1d.x
+    Valdes, September 24, 1986
+    1.  Changed subroutine name with a VOPS prefix to one with a FIT1D
+	prefix.
+
+pkg$images/imarith/t_imdivide.x
+pkg$images/doc/imdivide.hlp
+pkg$images/imdivide.par
+    Valdes, September 24, 1986
+    1.  Modified this ancient and obsolete task to remove redundant
+	subroutines now available in the VOPS library.
+    2.  The option to select action on zero divide was removed since
+	there was only one option.  Parameter file changed.
+    3.  Help page revised.
+
+pkg$images/geometry/t_blkrep.x +
+pkg$images/geometry/blkrp.gx +
+pkg$images/geometry/blkrep.x +
+pkg$images/doc/blkrep.hlp +
+pkg$images/doc/mkpkg
+pkg$images/images.cl
+pkg$images/images.men
+pkg$images/images.hd
+pkg$images/x_images.x
+    Valdes, September 24, 1986
+    1.  A new task called BLKREP for block replicating images has been added.
+	This task is a complement to BLKAVG and performs a function not
+	available in any other way.
+    2.  Help for BLKREP has been added.
+
+pkg$images/imarith/t_imarith.x
+pkg$images/imarith/imadiv.gx
+pkg$images/doc/imarith.hlp
+pkg$images/imarith.par
+    Valdes, September 24, 1986
+    1.  IMARITH has been modified to provide replacement of divisions
+	by zero with a constant parameter value.
+    2.  The documentation has been revised to include this change and to
+	clarify and emphasize areas of possible confusion.
+
+pkg$images/doc/magnify.hlp
+pkg$images/doc/blkavg.hlp
+    Valdes, September 18, 1986
+    1.  The MAGNIFY help document was expanded to clarify that images with axis
+        lengths of 1 cannot be magnified.  Also a discussion of the output
+	size of a magnified image.  This has been misunderstood often.
+    2.  Minor typo fix for BLKAVG.
+
+images$geometry/blkav.gx: Davis, September 7, 1986
+    1. The routine blkav$t was declared a function but called everywhere as
+    a procedure. Removed the function declaration.
+
+images$filters/med_sort.x: Davis, August 14, 1986
+    1. A bug in the sorting routine for MEDIAN and MODE in which the doop
+    loop increment was being set to zero has been fixed. This bug was
+    causing MEDIAN and MODE to fail on class 6 for certain sized windows.
+
+images$imfit/fit1d.x: Davis, July 24, 1986
+    1. A bug in the type=ratio option of fit1d was fixed. The iferr call
+    on the vector operator adivr was not trapping a divide by zero
+    condition. Changed adivr to adivzr.
+
+images$iminfo/listpixels.x: Davis, July 21, 1986
+    1. I changed a pargl to pargi for writing out the column number of the
+    pixels.
+
+images$iminfo/t_imstat.x: Davis, July 21, 1986
+    1. I changed a pargr to a pargd for the double precision quantitiies
+    sum(MIN) and sum(MAX).
+
+images$imfit/t_lineclean.x: Davis, July 14, 1986
+    1. Bug in the calling sequence for ic_clean fixed. The ic pointer
+    was not being passed to ic_clean causing access violation and/or
+    segmentation violation errors.
+
+images$imfit/fit1d.x, lineclean.x:  Valdes, July 3, 1986
+    1.  FIT1D and LINECLEAN modified to use new ICFIT package.
+
+From Valdes June 19, 1986
+
+1. The help page for IMSUM was modified to explicitly state what the
+median of an even number of images does.
+
+-----------------------------------------------------------------------------
+
+From Davis June 13, 1986
+
+1. A bug in CONVOLVE in which insufficient space was being allocated for
+long (> 161 elements) 1D kernels has been fixed. CONVOLVE was not
+allocating sufficent extra space.
+
+-----------------------------------------------------------------------------
+
+From Davis June 12, 1986
+
+1. I have changed the default value of parameter unmap in task FMEDIAN to
+yes to preserve the original data range.
+
+2. I have changed the value of parameter row_delimiter from \n to ;.
+
+-----------------------------------------------------------------------------
+
+From Davis May 12, 1986
+
+1. Changed the angle convention in GAUSS so that theta is the angle of the
+major axis with respect to the x axis measured counter-clockwise as specified
+in the help page instead of the negative of that angle.
+
+-----------------------------------------------------------------------------
+
+From Davis Apr 28, 1986
+
+1. Moved geomap.key to lib$scr and made redefined HELPFILE in geogmap.x
+appropriately.
+
+------------------------------------------------------------------------------
+
+images$imarith/imsum.gx:  Valdes Apr 25, 1986
+    1.  Fixed bug in generic code which called the real VOPS operator
+	regardless of the datatype.  This caused IMSUM to fail on short
+	images.
+
+From Davis Apr 17, 1986
+
+1. Changed constructs of the form boolean == false in the file imdelete.x
+to ! boolean.
+
+------------------------------------------------------------------------------
+
+images$imarith:  Valdes, April 8, 1986
+    1.  IMARITH has been modified to also operate on a list of specified
+	header parameters.  This is primarily used when adding images to
+	also added the exposure times.  A new parameter was added and the
+	help page modified.
+    2.  IMSUM has been modified to also operate on a list of specified
+	header parameters.  This is primarily used when summing images to
+	also sum the exposure times.  A new parameter was added and the
+	help page modified.
+
+------------------------------------------------------------------------------
+
+From Valdes Mar 24, 1986:
+
+1.  When modifying IMARITH to handle mixed dimensions the output image header
+was made a copy of the image with the higher dimension.  However, the default
+when the images were of the same dimension changed to be a copy of the second
+operand.  This has been changed back to being a copy of the first operand
+image.
+
+------------------------------------------------------------------------------
+
+From Davis Mar 21, 1986:
+
+1. A NULL pointer bug in the subroutine plfree inside IMSURFIT was causing
+segmentation violation errors. A null pointer test was added to plfree.
+
+------------------------------------------------------------------------------
+
+From Davis Mar 20, 1986:
+
+1. A bug involving in place operations in several image tasks has been  fixed.
+
+------------------------------------------------------------------------------
+
+From Davis Mar 19, 1986:
+
+1. IMSURFIT no longer permits the input image to be replaced by the output
+image.
+
+2. The tasks IMSHIFT, IMTRANSPOSE, SHIFTLINES, and GEOTRAN  have been modified
+to use the images tools xt_mkimtemp and xt_delimtemp for in place
+calculations.
+
+-------------------------------------------------------------------------------
+
+From Valdes Mar 13, 1986:
+
+1.  Bug dealing with type coercion in short datatype images in IMARITH and IMSUM
+which occurs on the SUN has been fixed.
+------
+From Valdes Mar 10, 1986:
+
+1.  IMSUM has been modified to work on any number of images.
+
+2.  Modified the help page
+------
+From Valdes Feb 25, 1986:
+
+There have been two changes to IMARITH:
+
+1.  A bug preventing use of image sections has been removed.
+
+2.  An improvement allowing use of images of different dimension.
+The algorithm is as follow:
+
+	a.  Check if both operands are images.  If not the output
+	image is a copy of the operand image.
+
+	b.  Check that the axes lengths are the same for the dimensions
+	in common.  For example a 3D and 2D image must have the same
+	number of columns and lines.
+
+	c.  Set the output image to be a copy of the image with the
+	higher dimension.
+
+	d.  Repeat the operation over the lower dimensions for each of
+	the higher dimensions.
+
+For example, consider subtracting a 2D image from a 3D image.  The output
+image will be 3D and the 2D image is subtracted from each band of the
+3D image.  This will work for any combination of dimensions.  Another
+example is dividing a 3D image by a 1D image.  Then each line of each
+plane and each band will be divided by the 1D image.  Likely applications
+will be subtracting biases and darks and dividing by response calibrations
+in stacked observations.
+
+3.  Modified the help page
+===========
+Release 2.2
+===========
+From Davis Mar 6, 1986:
+
+1. A serious bug had crept into GAUSS after I made some changes. For 2D
+images the sense of the sigma was reversed, i.e sigma = 2.0 was actually
+sigma = 0.5. This bug has now been fixed.
+
+---------------------------------------------------------------------------
+
+From Davis Jan 13, 1986:
+
+1. Listpixels will now print out complex pixel values correctly.
+
+---------------------------------------------------------------------------
+
+From Davis Dec 12, 1985:
+
+1. The directional gradient operator has been added to the images package.
+
+---------------------------------------------------------------------------
+
+From Valdes Dec 11, 1985:
+
+1.  IMARITH has been modified to first check if an operand is an existing
+file.  This allows purely numeric image names to be used.
+
+---------------------------------------------------------------------------
+
+From Davis Dec 11, 1985:
+
+1. A Laplacian (second derivatives) operator has been added to the images
+package.
+
+---------------------------------------------------------------------------
+
+From Davis Dec 10, 1985:
+
+1. The new convolution tasks  boxcar, gauss and convolve have been added
+to the images package. Convolve convolves an image with an arbitrary
+user supplied rectangular kernel. Gauss convolves an image with a 2D
+Gaussian of arbitrary size. Boxcar will smooth an image using a smoothing
+window of arbitrary size.
+
+2. The images package source code has been reorganized into the following
+subdirectories: 1) filters 2) geometry 3) imfit 4) imarith 4) iminfo and
+5) imutil 6) lib. Lib contains routines which may be of use to several IRAF
+tasks such as ranges. The imutil subdirectory contains tasks which modify
+images in some way such as hedit. The iminfo subdirectory contains code
+for displaying header and pixel values and other image characteristics
+such as the histogram. Image arithmetic and fitting routines are found
+in imarith and imfit respectively. Filters contains the convolution and
+median filtering routines and geometry contains the geometric distortion
+corrections routines.
+
+3. The documentation of the main images package has been brought into
+conformity with the new IRAF standards.
+
+4. Documentation for imdelete, imheader, imhistogram, listpixels and
+sections has been added to the help database.
+
+5. The parameter structure for imhistogram has been simplified. The
+redundant parameters sections and setranges have been removed.
+
+---------------------------------------------------------------------------
+
+
+From Valdes Nov 4, 1985:
+
+1.  IMCOPY modified so that the output image may be a directory.  Previously
+logical directories were not correctly identified.
+------
+
+From Davis Oct 21, 1985:
+
+1. A bug in the pixel rejection cycle of IMSURFIT was corrected. The routine
+make_ranges in ranges.x was not successfully converting a sorted list of
+rejected pixels into a list of ranges in all cases.
+
+2. Automatic zero divide error checking has been added to IMSURFIT.
+------
+From Valdes Oct 17, 1985:
+
+1.  Fit1d now allows averaging of image lines or columns when interactively
+setting the fitting parameters.  The syntax is "Fit line = 10 30"; i.e.
+blank separated line or column numbers.  A single number selects just one
+line or column.  Be aware however, that the actual fitting of the image
+is still done on each column or line individually.
+
+2.  The zero line in the interactive curve fitting graphs has been removed.
+This zero line interfered with fitting data near zero.
+------
+From Rooke Oct 10, 1985:
+
+1.  Blkaverage was changed to "blkavg" and modified to support any allowed
+number of dimensions.  It was also made faster in most cases, depending on 
+the blocking factors in each dimension.
+------
+From Valdes Oct 4, 1985:
+
+1.  Fit1d and lineclean modified to allow separate low and high rejection
+limits and rejection iterations.
+------
+From Davis Oct 3, 1985:
+
+1. Minmax was not calculating the minimum correctly for integer images.
+because the initial values were not being set correctly.
+------
+From Valdes Oct 1, 1985:
+
+1. Imheader was modified to print the image history.  Though the history
+mechanism is little used at the moment it should become an important part
+of any image.
+
+2.  Task revisions renamed to revs.
+------
+From Davis Sept 30, 1985:
+
+1. Two new tasks median and fmedian have been added to the images package.
+Fmedian is a fast median filtering algorithm for integer data which uses
+the histogram of the image to calculate the median at each window. Median
+is a slower but more general algorithm which performs the same task.
+------
+From Valdes August 26, 1985:
+
+1.  Blkaverage has been modified to include an new parameter called option.
+The current options are to average the blocks or sum the blocks.
+------
+From Valdes August 7, 1985
+
+1.  Fit1d and lineclean wer recompiled with the modified icfit package.
+The new package contains better labeling and graph documentation.
+
+2.  The two tasks now have parameters for setting the graphics device
+and reading cursor input from a file.
+______
+From: /u2/davis/ Tue 08:27:09 06-Aug-85
+Package: images
+Title: imshift bug
+ 
+Imshift was shifting incorrectly when an integral pixel shift in x and
+a fractional pixel shift in y was requested. The actual x shift was
+xshift + 1. The bug has been fixed and imshift will now work correctly for
+any combination of fractional and integral pixel shifts
+------
+From: /u2/davis/ Fri 18:14:12 02-Aug-85
+Package: images
+Title: new images task
+ 
+A new task GEOMAP has been added to the images package. GEOMAP calculates
+the spatial transformation required to map one image onto another.
+------
+From: /u2/davis/ Thu 16:47:49 01-Aug-85
+Package: images
+Title: new images tasks
+ 
+The tasks ROTATE, IMLINTRAN and GEODISTRAN have been added to the images
+package. ROTATE rotates and shifts an image. IMLINTRAN will rotate, rescale
+and shift an an image. GEODISTRAN corrects an image for geometric distortion.
+------
+From Valdes July 26, 1985:
+
+1.  The task revisions has been added to page revisions to the images
+package.  The intent is that each package will have a revisions task.
+Note that this means there may be multiple tasks named revisions loaded
+at one time.  Typing revisions alone will give the revisions for the
+current package.  To get the system revisions type system.revisions.
+
+2.  A new task called fit1d replaces linefit.  It is essentially the same
+as linefit except for an extra parameter "axis" which selects the axis along
+which the functions are to be fit.  Axis 1 is lines and axis 2 is columns.
+The advantages of this change are:
+
+	a.  Column fitting can now be done without transposing the image.
+	    This allows linefit to be used with image sections along
+	    both axes.
+	b.  For 1D images there is no prompt for the line number.
+.endhelp
diff --git a/pkg/images/imfit/doc/fit1d.hlp b/pkg/images/imfit/doc/fit1d.hlp
new file mode 100644
index 00000000..5b49f45f
--- /dev/null
+++ b/pkg/images/imfit/doc/fit1d.hlp
@@ -0,0 +1,177 @@
+.help fit1d Jul85 images.imfit
+.ih
+NAME
+fit1d -- fit a function to image lines
+.ih
+USAGE	
+.nf
+fit1d input output type
+.fi
+.ih
+PARAMETERS
+.ls input
+Images to be fit.  The images may contain image sections.  Only the region
+covered by the section will be modified in the output image.
+.le
+.ls output
+Output images to be created or modified.  The number of output images
+must match the number of input images.  If an output image does not exist
+it is first created and initialized to zero for fit types "fit" and
+"difference" and to one for fit type "ratio".
+.le
+.ls type
+Type of output.  The choices are:
+.ls fit      
+An image created from the function fits to the image lines.
+.le
+.ls difference
+The difference between the image and the fit (i.e. residuals).
+.le
+.ls ratio
+The ratio of the image and fit.
+.le
+.le
+.ls bpm = ""
+List of bad pixel masks.  This may be a null string to not use a
+bad pixel mask, a single mask that applies to all input images, or
+a matching list.  The value may also be !<keyword> to specify a keyword whose
+value is the mask to use.
+.le
+.ls axis = 1
+Axis along which the one dimensional fitting is done.  Axis 1 corresponds
+to fitting the image lines and axis 2 corresponds to fitting the columns.
+.le
+.ls interactive = yes
+If \fBinteractive\fR is set to yes, a plot of the fit is drawn and the
+cursor is available for interactively examining and adjusting the fit.
+.le
+.ls sample = "*"
+Lines or columns to be used in the fits.
+.le
+.ls naverage = 1
+Number of sample points to combined to create a fitting point.
+A positive value specifies an average and a negative value specifies
+a median.
+.le
+.ls function = spline3
+Function to be fit to the image lines or columns.  The functions are
+"legendre" (legendre polynomial), "chebyshev" (chebyshev polynomial),
+"spline1" (linear spline), and "spline3" (cubic spline).  The functions
+may be abbreviated.
+.le
+.ls order = 1
+The order of the polynomials or the number of spline pieces.
+.le
+.ls low_reject = 0., high_reject = 0.
+Rejection limits below and above the fit in units of the residual sigma.
+.le
+.ls niterate = 1
+Number of rejection iterations.
+.le
+.ls grow = 0.
+When a pixel is rejected, pixels within this distance of the rejected pixel
+are also rejected.
+.le
+.ls graphics = "stdgraph"
+Graphics output device for interactive graphics.
+.le
+.ls cursor = "stdgcur"
+Graphics cursor input.
+.le
+.ih
+DESCRIPTION
+A one dimensional function is fit to each line or column of the input images.
+The function may be a legendre polynomial, chebyshev polynomial,
+linear spline, or cubic spline of a given order or number of spline pieces.
+The output image is of pixel type real and is formed from the fitted
+function values, the difference or residuals of the fit (pixel value -
+fitted value), or the ratio between the pixel values and the fitted values.
+
+The output image may exist in which case a section in the input image is
+applied to the output image.  Thus, a section on the input image causes only
+that part of the output image to be changed.  If the output image does not
+exist it is first created with a size given by the full (without a section)
+input image and initialized to zero for fit and difference output types
+and one for ratio output types.
+
+A bad pixel mask may be specified to exclude data from the fitting.  Any
+non-zero value in the mask is excluded.   It appears in the interactive
+fitting in the same way as manually deleted points.  The mask is matched to
+the input image(s) as described by \fBpmmatch\fR.  The default is matching
+in physical coordinates.
+
+The points fit are determined by selecting a sample of lines or columns
+specified by the parameter \fIsample\fR and taking either the average or
+median of the number of points specified by the parameter \fInaverage\fR.
+The type of averaging is selected by the sign of the parameter and the number
+of points is selected by the absolute value of the parameter.
+The sample points are specified relative to any image sections.
+
+If \fIlow_reject\fR and/or \fIhigh_reject\fR are greater than zero the sigma
+of the residuals between the fitted points and the fitted function is computed
+and those points whose residuals are less than \fI-low_reject\fR * sigma
+and greater than \fIhigh_reject\fR * sigma are excluded from the fit.
+Points within a distance of \fIgrow\fR pixels of a rejected pixel are also
+excluded from the fit.  The function is then refit without the rejected points.
+This rejection procedure may be iterated a number of times given by the
+parameter \fIniterate\fR.
+
+The fitting parameters (\fIsample, naverage, function, order, low_reject,
+high_reject, niterate, grow\fR)
+may be adjusted interactively if the parameter \fIinteractive\fR is yes.
+Lines or columns from the image are selected to be fit with the \fBicfit\fR
+package.  A single column or line may be chosen or a blank-separated range
+may be averaged.  Note that the averaging applies only to the graphed
+data used to set the fitting parameters.  The actual image lines and columns
+are fit individually.  The interactive cursor mode commands for this package
+are described in a separate help entry under "icfit".  Line 1 is automatically
+selected for one dimensional images and any number of lines or columns may be
+selected for two dimensional images.  Note that the lines or columns are
+relative to the input image section; for example line 1 is the first line of
+the image section and not the first line of the image.  When an end-of-file or
+no line(s) or column(s) are given then the last selected fitting parameters
+are used on each line or column of the image.  This step is repeated for
+each image in the input list.
+.ih
+EXAMPLES
+1.  To create a smoothed version of an image by fitting the image lines:
+
+    cl> fit1d image fitimage fit
+
+If the interactive flag is set and the image is two dimensional then a prompt
+for an image line is printed:
+
+    image: Fit line = 100 200
+
+The selected lines are averaged, graphed, and the interactive options for
+setting and fitting the line are used.  Exiting with 'q' or return prompts for
+another line if the image is two dimensional.  When the fitting parameters
+are suitably set then respond with end-of-file or return to fit all the lines
+of the image and create the output image.
+
+2.  To subtract a linear function fit to columns 10 to 20 and 80 to 100 from
+columns 10 to 100 and to subtract another linear function fit to lines
+110 to 120 and 180 to 200 from columns 110 to 200:
+
+.nf
+    cl> fit1d image1[10:100,*] output diff axis=2 sample="1:11,71:91"
+    cl> fit1d image1[110:200,*] output diff axis=2 sample="1:11,71:91"
+.fi
+
+Pixels outside columns 10 to 100 and 110 to 200 are not affected.  Note that the
+sample points are specified relative to the image sections.  The script
+\fBbackground\fR is available in other packages for doing background
+subtractions.
+
+3.  To determine a small scale response image:
+
+    cl> fit1d image1 flat ratio
+
+The task \fBimred.generic.flat1d\fR is available for making flat field images
+by this method with the addition of an extra parameter to limit the data values
+for which the ratio is computed.
+.ih
+SEE ALSO
+imred.generic.background, imred.generic.flat1d
+xtools.icfit, lineclean, imsurfit
+.endhelp
diff --git a/pkg/images/imfit/doc/imsurfit.hlp b/pkg/images/imfit/doc/imsurfit.hlp
new file mode 100644
index 00000000..0aa46a97
--- /dev/null
+++ b/pkg/images/imfit/doc/imsurfit.hlp
@@ -0,0 +1,226 @@
+.help imsurfit Feb88 images.imfit
+.ih
+NAME
+imsurfit -- fit a surface function to an image
+.ih
+USAGE	
+imsurfit input, output, xorder, yorder
+.ih
+PARAMETERS
+.ls input
+List of images to be fit.
+.le
+.ls output
+Output image(s) of \fItype_output\fR.
+.le
+.ls xorder
+The order in x of the polynomials (1 = constant) or the number of polynomial
+pieces for the bicubic spline.
+.le
+.ls yorder
+The order in y of the polynomials (1 = constant) or the number of polynomial
+pieces for the bicubic spline.
+.le
+.ls cross_terms = yes
+Cross terms for the polynomials. For example, if \fIxorder\fR = 2 and
+\fIyorder\fR = 2
+then a function of the form z = a + b * x + c * y + d * x * y will be fit.
+.le
+.ls function = "leg"
+Functional for of surface to be fit to the image. The available functions
+(with minimum match abbreviation) are:
+.ls legendre
+.le
+.ls chebyshev
+.le
+.ls spline3
+.le
+.ls spline1
+.le
+.le
+.ls type_output = "fit"
+The type of output image.  The allowed types (with minimum match abbreviation)
+are:
+.ls clean
+The input image with deviant pixels in the good regions replaced by the
+fitted value.
+.le
+.ls fit  
+An image created from the surface fits to the image.
+.le
+.ls residual
+The difference of the input image and the fitted image.
+.le
+.ls response
+The ratio of the input image to the fitted image.
+All fitted (denominator) pixels below \fIdiv_min\fR are given a value of 1.
+.le
+.le
+.ls xmedian = 1, ymedian = 1
+The x and y dimensions of the box used for median processing.
+If \fIxmedian\fR > 1 or \fIymedian\fR is > 1,
+then the median is calculated for each box and used in the surface
+fit instead of the individual pixels.
+.le
+.ls median_percent = 50.
+If the number of pixels in the median box is less than \fImedian_percent\fR *
+\fIxmedian\fR * \fIymedian\fR the box will be omitted from the fit.
+.le
+.ls upper = 0., lower = 0.
+The number of sigma  limits for pixel rejection. If \fIupper\fR > 0. or
+\fIlower\fR > 0. and median processing is turned off,
+pixel rejection is enabled.
+.le
+.ls ngrow = 0
+The radius in pixels for region growing.
+Pixels within a distance of \fIngrow\fR pixels of
+a rejected pixel are also rejected.
+.le
+.ls niter = 0
+The maximum number of iterations in the rejection cycle.
+Rejection will be terminated if the number of rejected pixels is zero
+or the number of iterations equals \fIniter\fR.
+.le
+.ls regions = "all"
+The available options (with minimum match abbreviation) are:
+.ls all
+All points in the image are fit.
+.le
+.ls rows
+The fit is performed on the image rows specified by \fIrows\fR.
+.le
+.ls columns
+The fit is performed on the image columns specified by \fIcolumns\fR.
+.le
+.ls border
+The fit is performed on a border around the image whose width is specified
+by \fIborder\fR.
+.le
+.ls sections
+The fit is performed on image sections listed in the file specified
+by \fIsections\fR.
+.le
+.ls circle
+The fit is performed on a circular region whose parameters are specified by
+\fIcircle\fR.
+.le
+.ls invcircle
+The fit is performed on a region exterior to the circular region whose
+parameters are specified by \fIcircle\fR.
+.le
+.le
+.ls rows = "*"
+When \fIregion_type\fR = 'rows', the string parameter \fIrows\fR specifies
+the rows to be fit.
+.le
+.ls columns = "*"
+When \fIregion_type\fR = 'columns', the string parameter \fIcolumns\fR
+specifies the columns to be fit.
+.le
+.ls border = "50"
+When \fIregion_type\fR = 'border', the
+string parameter \fIborder\fR specifies the width of the border to be fit.
+.le
+.ls sections = ""
+When \fIregion_type\fR = 'sections', the
+string parameter \fIsections\fR is the name of the  file containing the list of
+image sections to be fit, where \fISections\fR may be the standard
+input (STDIN).
+The sections must be listed one per line in the following form: x1 x2 y1 y2.
+.le
+.ls circle = ""
+The string parameter \fIcircle\fR lists the parameter needed to specify
+the circle in the following format: xcenter ycenter radius. The three
+parameters must be integers.
+.le
+.ls div_min = INDEF
+When \fItype_output\fR = 'response' all divisions in which the fitted value
+is below \fIdiv_min\fR are set to the value 1.
+.le
+.ih
+DESCRIPTION
+A surface is fit to selected portions of the input image.
+The user may elect to fit the whole image, \fIregion_type\fR = 'all',
+selected rows, \fIregion_type\fR = 'rows', selected columns,
+\fIregion_type\fR = 'columns', a
+border around the image, \fIregion_type\fR = 'border' or image sections, 
+\fIregion_type\fR = 'sections'. If the sections  option is enabled the user
+must supply the name of the file containing a list of sections,
+\fIsections\fR = 'list', or enter them from the standard input. In either case
+the sections must be listed one per line in the following form: x1 x2 y1 y2.
+
+The parameter \fIsurface_type\fR may be a
+"legendre" polynomial, "chebyshev" polynomial,
+a cubic spline ("spline3") or a linear spline ("spline1").
+The order of the polynomials is selected in both x and y.
+Cross terms for the polynomial surfaces are optional.
+For the cubic spline the parameters \fIxorder\fR and \fIyorder\fR specify
+the number of polynomial pieces to be fit to the surface in
+each direction.
+
+The output image may be the fitted image, the difference between the
+input and the fitted image, the ratio of the input to the fitted image
+or the input image with deviant pixels in the fitted regions replaced
+with the fitted values, according to whether \fItype_output\fR =
+'fit', 'residual',
+'response' or 'clean'. If \fItype_output\fR = 'response' then pixels in the
+fitted image with values < \fIdiv_min\fR are replaced by 1.
+If \fItype_output\fR =
+'clean' then at least one of \fIupper\fR or \fIlower\fR must be > 0.
+
+Pixel rejection is enabled if median processing is turned off,
+\fIniter\fR > 0,
+and at least one of the parameters \fIupper\fR and \fIlower\fR is > 0.
+Region growing
+can be turned on by setting \fIngrow\fR > 0, in which case all pixels within
+a radius ngrow of a deviant pixel will be rejected.
+
+.ih
+EXAMPLES
+1. To create a smoothed version of an image:
+
+.nf
+	cl> imsurfit m74 m74smooth 5 10 function=spline3
+.fi
+
+2. To create a smoothed version of an image using median processing:
+
+.nf
+	cl> imsurfit m74 m74med 5 10 function=spline3 \
+	>>> xmed=5 ymed=5
+.fi
+
+3. To subtract a constant background from an image:
+
+.nf
+	cl> imsurfit abell30 abell30bck 1 1 function=leg \
+	>>> type=resid
+.fi
+
+4. To make a ratio image using signals above 1000 units:
+
+.nf
+	cl> imsurfit n7006 n7006ratio 20 20 function=spline3 \
+	>>> type=response div_min=1000
+.fi
+
+.ih
+TIMINGS
+Fitting and subtracting a constant from a 512 by 512 IRAF image requires
+~35 cpu seconds. Approximately 130 cpu seconds are required to fit a
+second degree polynomial in x and y (including cross-terms) to a
+100 pixel wide border around a 512 by
+512 IRAF image, and to subtract the fitted image from the input image.
+To produce a smooth 512 by 512 IRAF image using a 10 by 10 bicubic spline
+requires ~300 cpu seconds. Timings refer to a VAX 11/750 + fpa.
+
+.ih
+NOTES
+The surface fitting code uses the IRAF SURFIT math routines,
+which have been optimized for image fitting .
+The routines which fit selected portions
+of the image, perform pixel rejection and region growing, and create and
+maintain a list of rejected pixels utilize the ranges and pixlist packages
+of routines currently maintained in the images directory. These will be
+replaced by more general ranges and image masking routines in the future.
+.endhelp
diff --git a/pkg/images/imfit/doc/lineclean.hlp b/pkg/images/imfit/doc/lineclean.hlp
new file mode 100644
index 00000000..9ce2b95f
--- /dev/null
+++ b/pkg/images/imfit/doc/lineclean.hlp
@@ -0,0 +1,129 @@
+.help lineclean May85 images.imfit
+.ih
+NAME
+lineclean -- replace deviant pixels in image lines
+.ih
+USAGE	
+.nf
+lineclean input output
+.fi
+.ih
+PARAMETERS
+.ls input
+Input images to be cleaned.
+.le
+.ls output
+Output cleaned images.  The number of output images must be the same as the
+number of input images.
+.le
+.ls sample = "*"
+Columns to be used in fitting the cleaning function.
+.le
+.ls naverage = 1
+Number of sample points to combined to create a fitting point.
+A positive value specifies an average and a negative value specifies
+a median.
+.le
+.ls function = spline3
+Cleaning function to be fit to the image lines.  The functions are:
+.ls legendre
+Legendre polynomial of the specified order.
+.le
+.ls chebyshev
+Chebyshev polynomial of the specified order.
+.le
+.ls spline1
+Linear spline of the specified number of pieces.
+.le
+.ls spline3
+Cubic spline of the specified number of pieces.
+.le
+.le
+.ls order = 1
+The order of the polynomials or the number of spline pieces.
+.le
+.ls low_reject = 2.5, high_reject = 2.5
+Rejection limits below and above the fit in units of the residual sigma.
+.le
+.ls niterate = 1
+Number of rejection iterations.
+.le
+.ls grow = 1.
+When a pixel is rejected, pixels within this distance of the rejected pixel
+are also rejected.
+.le
+.ls graphics = "stdgraph"
+Graphics output device for interactive graphics.
+.le
+.ls cursor = "stdgcur"
+Graphics cursor input.
+.le
+.ih
+DESCRIPTION
+A one dimensional function is fit to each line of the input images.
+The function may be a legendre polynomial, chebyshev polynomial,
+linear spline, or cubic spline of a given order or number of spline pieces.
+If \fIlow_reject\fR and/or \fIhigh_reject\fR are greater than zero the sigma
+of the residuals between the fitted points and the fitted function is computed
+and those points whose residuals are less than \fI-low_reject\fR * sigma
+and greater than \fIhigh_reject\fR * sigma are excluded from the fit.
+Points within a distance of \fIgrow\fR pixels of a rejected pixel are also
+excluded from the fit.  The function is then refit without the rejected points.
+This rejection procedure may be iterated a number of times given by the
+parameter \fIniterate\fR.  Finally, the
+rejected points in the input image are replaced by the fitted values
+to create the output image lines.
+
+The output image may exist in which case a section in the input image is
+applied to the output image.  Thus, a section on the input image causes only
+that part of the output image to be cleaned.  If the output image does not
+exist it is first created by making a copy of the full (without a section)
+input image.
+
+The points fit are determined by selecting a sample of columns specified by
+the parameter \fIsample\fR and taking either the average or median of
+the number of points specified by the parameter \fInaverage\fR.
+The type of averaging is selected by the sign of the parameter and the number
+of points is selected by the absolute value of the parameter.
+The sample points are specified relative to any image section.
+
+The fitting parameters (\fIsample, naverage, function, order, low_reject,
+high_reject, niterate, grow\fR)
+may be adjusted interactively if the parameter \fIinteractive\fR is yes.
+Lines from the image are selected to be fit with the \fBicfit\fR package.
+For images of greater than two dimensions sets of numbers giving the
+2nd, 3rd, etc. coordinates are entered.
+The image lines are specified relative to any image section.
+When an end-of-file or no line is given then the last selected fitting
+parameters are used on each line of the image.  This step is repeated for
+each image in the input list.  The interactive options are described
+in the help information \fBicfit\fR.
+.ih
+EXAMPLES
+1. To clean pixels deviating by more than 2.5 sigma:
+
+	cl> lineclean image cleanimage
+
+If the interactive flag is set then a prompt for an image line is
+printed:
+
+	image: Fit line = 100
+
+For a one or two dimensional image the line number is entered (1 for a one
+dimensional image).  For a three dimensional image two numbers are entered.
+For example:
+
+	image: Fit line = 10 2
+
+for line 10 of the second image plane.
+
+The selected line is graphed and the interactive options for setting and
+fitting the line are used.  Data points marked with diamonds indicate
+points to be replaced by the fitted value.  Exiting with 'q' or return
+prompts for another line.  When the fitting parameters are suitably set
+then respond with end-of-file or return to fit all the lines of the image
+and create the output image.
+.ih
+SEE ALSO
+fit1d, xtools.icfit, imsurfit
+.endhelp
diff --git a/pkg/images/imfit/fit1d.par b/pkg/images/imfit/fit1d.par
new file mode 100644
index 00000000..f28100ba
--- /dev/null
+++ b/pkg/images/imfit/fit1d.par
@@ -0,0 +1,16 @@
+input,s,a,,,,Images to be fit
+output,s,a,,,,Output images
+bpm,s,h,"",,,Bad pixel mask(s)
+axis,i,h,1,1,2,Axis to be fit
+type,s,a,,,,"Type of output (fit, difference, ratio)"
+interactive,b,h,yes,,,Set fitting parameters interactively?
+sample,s,h,"*",,,Sample points to use in fit
+naverage,i,h,1,,,Number of points in sample averaging
+function,s,h,"spline3","spline3|legendre|chebyshev|spline1",,Fitting function
+order,i,h,1,1,,Order of fitting function
+low_reject,r,h,0.,0.,,Low rejection in sigma of fit
+high_reject,r,h,0.,0.,,High rejection in sigma of fit
+niterate,i,h,1,0,,Number of rejection iterations
+grow,r,h,0.,0.,,Rejection growing radius in pixels
+graphics,s,h,"stdgraph",,,Graphics output device
+cursor,*gcur,h,"",,,Graphics cursor input
diff --git a/pkg/images/imfit/imfit.cl b/pkg/images/imfit/imfit.cl
new file mode 100644
index 00000000..acde0687
--- /dev/null
+++ b/pkg/images/imfit/imfit.cl
@@ -0,0 +1,13 @@
+#{ IMFIT -- The Image Fitting Package.
+
+set	imfit		= "images$imfit/"
+
+package	imfit
+
+# Tasks.
+
+task	fit1d,
+	imsurfit,
+	lineclean	= "imfit$x_images.e"
+
+clbye()
diff --git a/pkg/images/imfit/imfit.hd b/pkg/images/imfit/imfit.hd
new file mode 100644
index 00000000..adcc4d44
--- /dev/null
+++ b/pkg/images/imfit/imfit.hd
@@ -0,0 +1,10 @@
+# Help directory for the IMFIT package
+
+$doc		= "images$imfit/doc/"
+$src		= "images$imfit/src/"
+
+fit1d		hlp=doc$fit1d.hlp,		src=src$fit1d.x
+imsurfit        hlp=doc$imsurfit.hlp,   	src=src$t_imsurfit.x
+lineclean       hlp=doc$lineclean.hlp,  	src=src$t_lineclean.x
+revisions	sys=Revisions
+
diff --git a/pkg/images/imfit/imfit.men b/pkg/images/imfit/imfit.men
new file mode 100644
index 00000000..793b1c5e
--- /dev/null
+++ b/pkg/images/imfit/imfit.men
@@ -0,0 +1,3 @@
+          fit1d - Fit a function to image lines or columns
+       imsurfit - Fit a surface to a 2-D image
+      lineclean - Replace deviant pixels in image lines
diff --git a/pkg/images/imfit/imfit.par b/pkg/images/imfit/imfit.par
new file mode 100644
index 00000000..cef3f3ff
--- /dev/null
+++ b/pkg/images/imfit/imfit.par
@@ -0,0 +1 @@
+version,s,h,"Jan97"
diff --git a/pkg/images/imfit/imsurfit.par b/pkg/images/imfit/imsurfit.par
new file mode 100644
index 00000000..1f21ffe5
--- /dev/null
+++ b/pkg/images/imfit/imsurfit.par
@@ -0,0 +1,24 @@
+# IMSURFIT
+
+input,f,a,,,,Input images to be fit
+output,f,a,,,,Output images
+xorder,i,a,2,1,,Order of function in x
+yorder,i,a,2,1,,Order of function in y
+type_output,s,h,'fit',,,'Type of output (fit,residual,response,clean)'
+function,s,h,'leg',,,'Function to be fit (legendre,chebyshev,spline3)'
+cross_terms,b,h,y,,,Include cross-terms for polynomials?
+xmedian,i,h,1,1,,X length of median box
+ymedian,i,h,1,1,,Y length of median box
+median_percent,r,h,50.,,,Minimum fraction of pixels in median box
+lower,r,h,0.0,0.0,,Lower limit for residuals
+upper,r,h,0.0,0.0,,Upper limit for residuals
+ngrow,i,h,0,0,,Radius of region growing circle
+niter,i,h,0,0,,Maximum number of rejection cycles
+regions,s,h,'all',,, 'Good regions (all,rows,columns,border,sections,circle,invcircle)'
+rows,s,h,'*',,,Rows to be fit
+columns,s,h,'*',,,Columns to be fit
+border,s,h,'50',,,Width of border to be fit
+sections,s,h,,,,File name for sections list
+circle,s,h,,,,Circle specifications
+div_min,r,h,INDEF,,,Division minimum for response output
+mode,s,h,'ql'
diff --git a/pkg/images/imfit/lineclean.par b/pkg/images/imfit/lineclean.par
new file mode 100644
index 00000000..5942f03f
--- /dev/null
+++ b/pkg/images/imfit/lineclean.par
@@ -0,0 +1,13 @@
+input,s,a,,,,Images to be cleaned
+output,s,a,,,,Output images
+interactive,b,h,yes,,,Set fitting parameters interactively?
+sample,s,h,"*",,,Sample points to use in fit
+naverage,i,h,1,,,Number of points in sample averaging
+function,s,h,"spline3","spline3|legendre|chebyshev|spline1",,Fitting function
+order,i,h,1,1,,Order of fitting function
+low_reject,r,h,2.5,0.,,Low rejection in sigma of fit
+high_reject,r,h,2.5,0.,,High rejection in sigma of fit
+niterate,i,h,1,0,,Number of rejection iterations
+grow,r,h,1.,0.,,Rejection growing radius in pixels
+graphics,s,h,"stdgraph",,,Graphics output device
+cursor,*gcur,h,"",,,Graphics cursor input
diff --git a/pkg/images/imfit/mkpkg b/pkg/images/imfit/mkpkg
new file mode 100644
index 00000000..a74a627e
--- /dev/null
+++ b/pkg/images/imfit/mkpkg
@@ -0,0 +1,5 @@
+# MKPKG for the IMFIT Package
+
+libpkg.a:
+	@src
+	;
diff --git a/pkg/images/imfit/src/fit1d.x b/pkg/images/imfit/src/fit1d.x
new file mode 100644
index 00000000..84ffddf1
--- /dev/null
+++ b/pkg/images/imfit/src/fit1d.x
@@ -0,0 +1,597 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include <pkg/gtools.h>
+include	<error.h>
+
+define	MAXBUF	(512*100)		# Maximum number of pixels per block
+
+
+# FIT1D -- Fit a function to image lines or columns and output an image
+# consisting of the fit, the difference, or the ratio.  The fitting parameters
+# may be set interactively using the icfit package.
+
+procedure t_fit1d ()
+
+int	listin				# Input image list
+int	listout				# Output image list
+int	listbpm				# Bad pixel mask list
+bool	interactive			# Interactive?
+
+char	sample[SZ_LINE]			# Sample ranges
+int	naverage			# Sample averaging size
+char	function[SZ_LINE]		# Curve fitting function
+int	order				# Order of curve fitting function
+real	low_reject, high_reject		# Rejection thresholds
+int	niterate			# Number of rejection iterations
+real	grow				# Rejection growing radius
+
+int	axis				# Image axis to fit
+int	ntype				# Type of output
+char	input[SZ_LINE]			# Input image
+char	output[SZ_FNAME]		# Output image
+char	bpm[SZ_FNAME]			# Bad pixel mask
+pointer	in, out, bp			# IMIO pointers
+pointer	ic				# ICFIT pointer
+pointer	gt				# GTOOLS pointer
+
+bool	same, clgetb()
+int	imtopen(), imtgetim(), imtlen(), strdic(), gt_init()
+int	clgeti()
+real	clgetr()
+
+begin
+	# Get input and output lists and check that the number of images
+	# are the same.
+
+	call clgstr ("input", input, SZ_LINE)
+	listin = imtopen (input)
+	call clgstr ("output", input, SZ_LINE)
+	listout = imtopen (input)
+	if (imtlen (listin) != imtlen (listout)) {
+	    call imtclose (listin)
+	    call imtclose (listout)
+	    call error (0, "Input and output image lists do not match")
+	}
+	call clgstr ("bpm", input, SZ_LINE)
+	listbpm = imtopen (input)
+	if (imtlen (listbpm) > 1 && imtlen (listin) != imtlen (listbpm)) {
+	    call imtclose (listin)
+	    call imtclose (listout)
+	    call imtclose (listbpm)
+	    call error (0, "Input and mask lists do not match")
+	}
+
+	# Get task parameters.
+
+	axis = clgeti ("axis")
+	call clgstr ("type", input, SZ_LINE)
+	call clgstr ("sample", sample, SZ_LINE)
+	naverage = clgeti ("naverage")
+	call clgstr ("function", function, SZ_LINE)
+	order = clgeti ("order")
+	low_reject = clgetr ("low_reject")
+	high_reject = clgetr ("high_reject")
+	niterate = clgeti ("niterate")
+	grow = clgetr ("grow")
+	interactive = clgetb ("interactive")
+
+ 	# Decode the output type and initialize the curve fitting package.
+
+	ntype = strdic (input, input, SZ_LINE, "|fit|difference|ratio|")
+	if (ntype == 0)
+	    call error (0, "Unknown output type")
+
+	# Set the ICFIT pointer structure.
+	call ic_open (ic)
+	call ic_pstr (ic, "sample", sample)
+	call ic_puti (ic, "naverage", naverage)
+	call ic_pstr (ic, "function", function)
+	call ic_puti (ic, "order", order)
+	call ic_putr (ic, "low", low_reject)
+	call ic_putr (ic, "high", high_reject)
+	call ic_puti (ic, "niterate", niterate)
+	call ic_putr (ic, "grow", grow)
+	call ic_pstr (ic, "ylabel", "")
+
+	gt = gt_init()
+	call gt_sets (gt, GTTYPE, "line")
+
+	# Fit the lines in each input image.
+
+	bpm[1] = EOS
+	while ((imtgetim (listin, input, SZ_LINE) != EOF) &&
+	    (imtgetim (listout, output, SZ_FNAME) != EOF)) {
+	    if (imtgetim (listbpm, bpm, SZ_FNAME) == EOF)
+		;
+
+	    iferr (call f1d_immap (input,output,bpm,ntype,in,out,bp,same)) {
+		call erract (EA_WARN)
+		next
+	    }
+	    call f1d_fit1d (in,out,bp,ic,gt,input,axis,ntype,interactive)
+	    call imunmap (in)
+	    if (!same)
+		call imunmap (out)
+	    if (bp != NULL)
+		call yt_pmunmap (bp)
+	}
+
+	call ic_closer (ic)
+	call gt_free (gt)
+	call imtclose (listin)
+	call imtclose (listout)
+end
+
+
+# F1D_FIT1D -- Given the image descriptor determine the fitting function
+# for each line or column and create an output image.  If the interactive flag
+# is set then set the fitting parameters interactively.
+
+procedure f1d_fit1d (in, out, bp, ic, gt, title, axis, ntype, interactive)
+
+pointer	in				# IMIO pointer for input image
+pointer	out				# IMIO pointer for output image
+pointer	bp				# IMIO pointer for bad pixel mask
+pointer	ic				# ICFIT pointer
+pointer	gt				# GTOOLS pointer
+char	title[ARB]			# Title
+int	axis				# Image axis to fit
+int	ntype				# Type of output
+bool	interactive			# Interactive?
+
+char	graphics[SZ_FNAME]		# Graphics device
+int	i, nx, new
+real	div
+pointer	cv, gp, sp, x, wts, indata, outdata
+
+int	f1d_getline(), f1d_getdata(), strlen()
+real	cveval()
+pointer	gopen()
+
+begin
+	# Error check.
+
+	if (IM_NDIM (in) > 2)
+	    call error (0, "Image dimensions > 2 are not implemented")
+	if (axis > IM_NDIM (in))
+	    call error (0, "Axis exceeds image dimension")
+
+	# Allocate memory for curve fitting.
+
+	nx = IM_LEN (in, axis)
+	call smark (sp)
+	call salloc (x, nx, TY_REAL)
+
+	do i = 1, nx
+	    Memr[x+i-1] = i
+
+	call ic_putr (ic, "xmin", Memr[x])
+	call ic_putr (ic, "xmax", Memr[x+nx-1])
+
+	# If the interactive flag is set then use icg_fit to set the
+	# fitting parameters.  Get_fitline returns EOF when the user
+	# is done.  The weights are reset since the user may delete
+	# points.
+
+	if (interactive) {
+	    call clgstr ("graphics", graphics, SZ_FNAME)
+	    gp = gopen (graphics, NEW_FILE, STDGRAPH)
+
+	    i = strlen (title)
+	    indata = NULL
+	    while (f1d_getline (ic,gt,in,bp,axis,title,indata,wts) != EOF) {
+		title[i + 1] = EOS
+	        call icg_fit (ic, gp, "cursor", gt, cv, Memr[x], Memr[indata],
+		    Memr[wts], nx)
+	    }
+	    call mfree (indata, TY_REAL)
+	    call mfree (wts, TY_REAL)
+	    call gclose (gp)
+	}
+
+	# Loop through the input image and create an output image.
+
+	new = YES
+
+	while (f1d_getdata (in,out,bp,axis,MAXBUF,indata,outdata,wts) != EOF) {
+
+	    call ic_fit (ic, cv, Memr[x], Memr[indata], Memr[wts],
+		nx, new, YES, new, new)
+	    new = NO
+
+	    # Be careful because the indata and outdata buffers may be the same.
+	    switch (ntype) {
+	    case 1:
+	        call cvvector (cv, Memr[x], Memr[outdata], nx)
+	    case 2:
+		do i = 0, nx-1
+		    Memr[outdata+i] = Memr[indata+i] - cveval (cv, Memr[x+i])
+	    case 3:
+		do i = 0, nx-1 {
+		    div = cveval (cv, Memr[x+i])
+		    if (abs (div) < 1E-20)
+			div = 1
+		    Memr[outdata+i] = Memr[indata+i] / div
+		}
+	    }
+	}
+
+	call cvfree (cv)
+	call mfree (wts, TY_REAL)
+	call sfree (sp)
+end
+
+
+# F1D_IMMAP -- Map images for fit1d.
+
+procedure f1d_immap (input, output, bpm, ntype, in, out, bp, same)
+
+char	input[ARB]		# Input image
+char	output[ARB]		# Output image
+char	bpm[ARB]		# Bad pixel mask
+int	ntype			# Type of fit1d output
+pointer	in			# Input IMIO pointer
+pointer	out			# Output IMIO pointer
+pointer	bp			# Mask IMIO pointer
+bool	same			# Same image?
+
+int	i
+pointer	sp, iroot, isect, oroot, osect, line, data
+
+bool	streq()
+int	imaccess(), impnlr()
+pointer	immap(), yt_pmmap()
+errchk	immap, yt_pmmap
+
+begin
+	# Get the root name and section of the input image.
+
+	call smark (sp)
+	call salloc (iroot, SZ_FNAME, TY_CHAR)
+	call salloc (isect, SZ_FNAME, TY_CHAR)
+	call salloc (oroot, SZ_FNAME, TY_CHAR)
+	call salloc (osect, SZ_FNAME, TY_CHAR)
+
+	call imgimage (input, Memc[iroot], SZ_FNAME)
+	call imgsection (input, Memc[isect], SZ_FNAME)
+	call imgimage (output, Memc[oroot], SZ_FNAME)
+	call imgsection (output, Memc[osect], SZ_FNAME)
+	same = streq (Memc[iroot], Memc[oroot])
+
+	# If the output image is not accessible then create it as a new copy
+	# of the full input image and initialize according to ntype.
+
+	if (imaccess (output, READ_WRITE) == NO) {
+	    in = immap (Memc[iroot], READ_ONLY, 0)
+	    out = immap (Memc[oroot], NEW_COPY, in)
+	    IM_PIXTYPE(out) = TY_REAL
+
+	    call salloc (line, IM_MAXDIM, TY_LONG)
+	    call amovkl (long (1), Meml[line], IM_MAXDIM)
+
+	    switch (ntype) {
+	    case 1, 2:
+	        while (impnlr (out, data, Meml[line]) != EOF)
+	            call aclrr (Memr[data], IM_LEN(out, 1))
+	    case 3:
+	        while (impnlr (out, data, Meml[line]) != EOF)
+	            call amovkr (1., Memr[data], IM_LEN(out, 1))
+	    }
+
+	    call imunmap (in)
+	    call imunmap (out)
+	}
+
+	# Map the images.  If the output image has a section
+	# then use it.  If the input image has a section and the output image
+	# does not then add the image section to the output image.  Finally
+	# check the input and output images have the same size.
+
+	in = immap (input, READ_ONLY, 0)
+
+	if (Memc[isect] != EOS && Memc[osect] == EOS) {
+	    call sprintf (Memc[osect], SZ_FNAME, "%s%s")
+	        call pargstr (Memc[oroot])
+	        call pargstr (Memc[isect])
+	} else
+	    call strcpy (output, Memc[osect], SZ_FNAME)
+
+	if (streq (input, Memc[osect])) {
+	    call imunmap (in)
+	    in = immap (input, READ_WRITE, 0)
+	    out = in
+	} else
+	    out = immap (Memc[osect], READ_WRITE, 0)
+
+	do i = 1, IM_NDIM(in)
+	    if (IM_LEN(in, i) != IM_LEN(out, i)) {
+		call imunmap (in)
+		if (!same)
+		    call imunmap (out)
+		call sfree (sp)
+		call error (0, "Input and output images have different sizes")
+	    }
+
+	bp = yt_pmmap (bpm, in, Memc[iroot], SZ_FNAME)
+
+	call sfree (sp)
+end
+
+
+# F1D_GETDATA -- Get a line of image data.
+
+int procedure f1d_getdata (in, out, bp, axis, maxbuf, indata, outdata, wts)
+
+pointer	in			# Input IMIO pointer
+pointer	out			# Output IMIO pointer
+pointer	bp			# Bad pixel mask IMIO pointer
+int	axis			# Image axis
+int	maxbuf			# Maximum buffer size for column axis
+pointer	indata			# Input data pointer
+pointer	outdata			# Output data pointer
+pointer	wts			# Weights pointer
+
+int	i, index, last_index, col1, col2, nc, ncols, nlines, ncols_block
+pointer	inbuf, outbuf, wtbuf, ptr
+pointer	imgl1r(), imgl1s(), impl1r(), imgl2r(), imgl2s(), impl2r()
+pointer	imgs2r(), imgs2s(), imps2r()
+data	index/0/
+
+begin
+	# Increment to the next image vector.
+
+	index = index + 1
+
+	# Initialize for the first vector.
+
+	if (index == 1) {
+	    ncols = IM_LEN (in, 1)
+	    if (IM_NDIM (in) == 1)
+		nlines = 1
+	    else
+		nlines = IM_LEN (in, 2)
+
+	    switch (axis) {
+	    case 1:
+		last_index = nlines
+
+	        call malloc (wts, ncols, TY_REAL)
+	    case 2:
+		last_index = ncols
+	        ncols_block = max (1, min (ncols, maxbuf / nlines))
+		col2 = 0
+
+	        call malloc (indata, nlines, TY_REAL)
+	        call malloc (outdata, nlines, TY_REAL)
+	        call malloc (wts, nlines, TY_REAL)
+	    }
+	}
+
+	# Finish up if the last vector has been done.
+
+	if (index > last_index) {
+	    switch (axis) {
+	    case 1:
+	        call mfree (wts, TY_REAL)
+	    case 2:
+	        ptr = outbuf + index - 1 - col1
+	        do i = 1, nlines {
+		    Memr[ptr] = Memr[outdata+i-1]
+		    ptr = ptr + nc
+	        }
+
+	        call mfree (indata, TY_REAL)
+	        call mfree (outdata, TY_REAL)
+	        call mfree (wts, TY_REAL)
+	    }
+
+	    index = 0
+	    return (EOF)
+	}
+
+	# Get the next image vector.
+
+	switch (axis) {
+	case 1:
+	    ncols = IM_LEN(in,1)
+	    if (IM_NDIM (in) == 1) {
+		indata = imgl1r (in)
+		outdata = impl1r (out)
+		if (bp == NULL)
+		    call amovkr (1., Memr[wts], ncols)
+		else {
+		    wtbuf = imgl1s (bp)
+		    do i = 0, ncols-1 {
+		        if (Mems[wtbuf+i] == 0)
+			    Memr[wts+i] = 1.
+			else
+			    Memr[wts+i] = 0.
+		    }
+		}
+	    } else {
+		indata = imgl2r (in, index)
+		outdata = impl2r (out, index)
+		if (bp == NULL)
+		    call amovkr (1., Memr[wts], ncols)
+		else {
+		    wtbuf = imgl2s (bp, index)
+		    do i = 0, ncols-1 {
+		        if (Mems[wtbuf+i] == 0)
+			    Memr[wts+i] = 1.
+			else
+			    Memr[wts+i] = 0.
+		    }
+		}
+	    }
+	case 2:
+	    if (index > 1) {
+		ptr = outbuf + index - 1 - col1
+		do i = 1, nlines {
+		    Memr[ptr] = Memr[outdata+i-1]
+		    ptr = ptr + nc
+		}
+	    }
+
+	    if (index > col2) {
+		col1 = col2 + 1
+		col2 = min (ncols, col1 + ncols_block - 1)
+		inbuf = imgs2r (in, col1, col2, 1, nlines)
+		outbuf = imps2r (out, col1, col2, 1, nlines)
+		if (bp != NULL)
+		    wtbuf = imgs2s (bp, col1, col2, 1, nlines)
+		nc = col2 - col1 + 1
+	    }
+
+	    ptr = inbuf + index - col1
+	    do i = 0, nlines-1 {
+		Memr[indata+i] = Memr[ptr]
+		ptr = ptr + nc
+	    }
+
+	    if (bp == NULL)
+		call amovkr (1., Memr[wts], nlines)
+	    else {
+		ptr = wtbuf + index - col1
+		do i = 0, nlines-1 {
+		    if (Mems[ptr] == 0)
+			Memr[wts+i] = 1.
+		    else
+			Memr[wts+i] = 0.
+		    ptr = ptr + nc
+		}
+	    }
+	}
+
+	return (index)
+end
+
+
+# F1D_GETLINE -- Get image data to be fit interactively.  Return EOF
+# when the user enters EOF or CR.  Default is 1 and the out of bounds
+# requests are silently limited to the nearest in edge.
+
+int procedure f1d_getline (ic, gt, im, bp, axis, title, data, wts)
+
+pointer	ic			# ICFIT pointer
+pointer	gt			# GTOOLS pointer
+pointer	im			# IMIO pointer input image
+pointer	bp			# IMIO pointer for bad pixel mask
+int	axis			# Image axis
+char	title[ARB]		# Title
+pointer	data			# Image data
+pointer	wts			# Weights
+
+pointer	x, wtbuf
+char	line[SZ_LINE]
+int	i, j, stat, imlen
+int	getline(), nscan()
+pointer	imgl1r(), imgl1s()
+data	stat/EOF/
+
+begin
+	# If the image is one dimensional do not prompt.
+
+	if (IM_NDIM (im) == 1) {
+	    if (stat == EOF) {
+		call sprintf (title, SZ_LINE, "%s\n%s")
+	    	    call pargstr (title)
+	    	    call pargstr (IM_TITLE(im))
+		call gt_sets (gt, GTTITLE, title)
+		call mfree (data, TY_REAL)
+		imlen = IM_LEN(im,1)
+		call malloc (data, imlen, TY_REAL)
+		call amovr (Memr[imgl1r(im)], Memr[data], imlen)
+		call malloc (wts, imlen, TY_REAL)
+		if (bp == NULL)
+		    call amovkr (1., Memr[wts], imlen)
+		else {
+		    wtbuf = imgl1s (bp)
+		    do i = 0, imlen-1 {
+		        if (Mems[wtbuf+i] == 0)
+			    Memr[wts+i] = 1.
+			else
+			    Memr[wts+i] = 0.
+		    }
+		}
+		stat = OK
+	    } else
+		stat = EOF
+
+	    return (stat)
+	}
+
+	# If the image is two dimensional prompt for the line or column.
+
+	switch (axis) {
+	case 1:
+	    imlen = IM_LEN (im, 2)
+	    call sprintf (title, SZ_LINE, "%s: Fit line =")
+	        call pargstr (title)
+	case 2:
+	    imlen = IM_LEN (im, 1)
+	    call sprintf (title, SZ_LINE, "%s: Fit column =")
+	        call pargstr (title)
+	}
+
+	call printf ("%s ")
+	    call pargstr (title)
+	call flush (STDOUT)
+
+	if (getline(STDIN, line) == EOF)
+	    return (EOF)
+
+	call sscan (line)
+	call gargi (i)
+	call gargi (j)
+
+	switch (nscan()) {
+	case 0:
+	    stat = EOF
+	    return (stat)
+	case 1:
+	    i = max (1, min (imlen, i))
+	    j = i
+	case 2:
+	    i = max (1, min (imlen, i))
+	    j = max (1, min (imlen, j))
+	}
+
+	call sprintf (title, SZ_LINE, "%s %d - %d\n%s")
+	    call pargstr (title)
+	    call pargi (i)
+	    call pargi (j)
+	    call pargstr (IM_TITLE(im))
+
+	call gt_sets (gt, GTTITLE, title)
+
+	call malloc (wts, imlen, TY_REAL)
+	switch (axis) {
+	case 1:
+	    call ic_pstr (ic, "xlabel", "Column")
+	    call xt_21imavg (im, axis, 1, IM_LEN(im,1), i, j, x, data, imlen)
+	    if (bp != NULL)
+		call xt_21imsum (bp, axis, 1, IM_LEN(im,1), i, j, x, wts, imlen)
+	case 2:
+	    call ic_pstr (ic, "xlabel", "Line")
+	    call xt_21imavg (im, axis, i, j, 1, IM_LEN(im,2), x, data, imlen)
+	    if (bp != NULL)
+		call xt_21imsum (bp, axis, i, j, 1, IM_LEN(im,2), x, wts, imlen)
+	}
+	if (bp == NULL) {
+	    call mfree (wts, TY_REAL)
+	    call malloc (wts, imlen, TY_REAL)
+	    call amovkr (1., Memr[wts], imlen)
+	} else {
+	    do i = 0, imlen-1 {
+		if (Memr[wts+i] == 0.)
+		    Memr[wts+i] = 1.
+		else
+		    Memr[wts+i] = 0.
+	    }
+	}
+	call mfree (x, TY_REAL)
+
+	stat = OK
+	return (stat)
+end
diff --git a/pkg/images/imfit/src/imsurfit.h b/pkg/images/imfit/src/imsurfit.h
new file mode 100644
index 00000000..84c077ec
--- /dev/null
+++ b/pkg/images/imfit/src/imsurfit.h
@@ -0,0 +1,40 @@
+# Header file for IMSURFIT
+
+define	LEN_IMSFSTRUCT	20
+
+# surface parameters
+define	SURFACE_TYPE	Memi[$1]
+define	XORDER		Memi[$1+1]
+define	YORDER		Memi[$1+2]
+define	CROSS_TERMS	Memi[$1+3]
+define	TYPE_OUTPUT	Memi[$1+4]
+
+# median processing parameters
+define	MEDIAN		Memi[$1+5]
+define	XMEDIAN		Memi[$1+6]
+define	YMEDIAN		Memi[$1+7]
+define	MEDIAN_PERCENT	Memr[P2R($1+8)]
+
+# pixel rejection parameters
+define	REJECT		Memi[$1+9]
+define	NGROW		Memi[$1+10]
+define	NITER		Memi[$1+11]
+define	LOWER		Memr[P2R($1+12)]
+define	UPPER		Memr[P2R($1+13)]
+
+define	DIV_MIN		Memr[P2R($1+14)]
+
+# definitions for type_output
+define	FIT		1
+define	CLEAN		2
+define	RESID		3
+define	RESP		4
+
+# definitions for good regions
+define	ALL		1
+define	COLUMNS		2
+define	ROWS		3
+define	BORDER		4
+define	SECTIONS	5
+define	CIRCLE		6
+define	INVCIRCLE	7
diff --git a/pkg/images/imfit/src/imsurfit.x b/pkg/images/imfit/src/imsurfit.x
new file mode 100644
index 00000000..9f655f52
--- /dev/null
+++ b/pkg/images/imfit/src/imsurfit.x
@@ -0,0 +1,1172 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <mach.h>
+include <imhdr.h>
+include <math/surfit.h>
+include "imsurfit.h"
+
+# IMSURFIT -- Procedure to fit a surface to a single image including
+# optional pixel rejection.
+
+procedure imsurfit (imin, imout, imfit, gl)
+
+pointer	imin	# pointer to the input image
+pointer	imout	# pointer to the output image
+pointer	imfit	# pointer to the imsurfit parameters
+pointer	gl	# pointer to the good regions list
+
+pointer	sf, rl
+errchk	isfree, prl_free
+errchk	all_pixels, good_pixels, good_median, all_medians, do_reject
+errchk	set_outimage
+
+begin
+	sf = NULL
+	rl = NULL
+
+	# Accumulate and solve the surface.
+	if (gl == NULL) {
+	    if (MEDIAN(imfit) == NO)
+	        call all_pixels (imin, imfit, sf)
+	    else
+		call all_medians (imin, imfit, sf)
+	} else {
+	    if (MEDIAN(imfit) == NO)
+	        call good_pixels (imin, imfit, gl, sf)
+	    else
+		call good_medians (imin, imfit, gl, sf)
+	}
+
+	# Perform the reject cycle.
+	if (REJECT(imfit) == YES || TYPE_OUTPUT(imfit) == CLEAN)
+	    call do_reject (imin, imfit, gl, sf, rl)
+
+	# Evaluate surface for appropriate output type.
+	call set_outimage (imin, imout, imfit, sf, rl)
+
+	# Cleanup.
+	call prl_free (rl)
+	call isfree (sf)
+
+	rl = NULL
+	sf = NULL
+end
+
+
+# ALL_PIXELS -- Accumulate surface when there are no bad regions
+# and no median processing.
+
+procedure all_pixels (im, imfit, sf)
+
+pointer	im	# pointer to the input image
+pointer	imfit	# pointer to the imsurfit structure
+pointer	sf	# pointer to the surface descriptor
+
+int	i, lp, ncols, nlines, ier
+long	v[IM_MAXDIM]
+pointer	sp, cols, lines, wgt, lbuf
+int	imgnlr()
+errchk	smark, salloc, sfree, imgnlr
+errchk isinit, islfit, islrefit, issolve
+
+begin
+	ncols = IM_LEN(im, 1)
+	nlines = IM_LEN(im,2)
+
+	# Initialize the surface fit.
+	call isinit (sf, SURFACE_TYPE(imfit), XORDER(imfit), YORDER(imfit),
+	    CROSS_TERMS(imfit), ncols, nlines) 
+
+	# Allocate working space for fitting.
+	call smark (sp)
+	call salloc (cols, ncols, TY_INT)
+	call salloc (lines, nlines, TY_INT)
+	call salloc (wgt, ncols, TY_REAL)
+
+	# Initialize the x and weight buffers.
+	do i = 1, ncols
+	    Memi[cols - 1 + i] = i
+	 call amovkr (1.0, Memr[wgt], ncols)
+
+	# Loop over image lines.
+	lp = 0
+	call amovkl (long (1), v, IM_MAXDIM)
+	do i = 1, nlines {
+
+	    # Read in the image line.
+	    if (imgnlr (im, lbuf, v) == EOF)
+		call error (0, "Error reading image")
+
+	    # Fit each image line.
+	    if (i == 1)
+	        call islfit (sf, Memi[cols], i, Memr[lbuf], Memr[wgt],
+		ncols, SF_USER, ier)
+	    else
+		call islrefit (sf, Memi[cols], i, Memr[lbuf], Memr[wgt])
+
+	    # Handle fitting errors.
+	    switch (ier) {
+	    case NO_DEG_FREEDOM:
+		call eprintf ("Warning: Too few columns to fit line: %d\n")
+		    call pargi (i)
+	    case SINGULAR:
+		call eprintf ("Warning: Solution singular for line: %d\n")
+		    call pargi (i)
+		Memi[lines + lp] = i
+		lp = lp + 1
+	    default:
+		Memi[lines + lp] = i
+		lp = lp + 1
+	    }
+
+	}
+
+	# Solve the surface.
+	call issolve (sf, Memi[lines], lp, ier)
+
+	# Handle fitting errors.
+	switch (ier) {
+	case NO_DEG_FREEDOM:
+	    call error (0, "ALL_PIXELS: Cannot fit surface.")
+	case SINGULAR:
+	    call eprintf ("Warning: Solution singular for surface.\n")
+	default:
+	    # everything OK
+	}
+
+	# Free space.
+	call sfree (sp)
+end
+
+
+# GOOD_PIXELS -- Get surface when good regions are defined and median
+# processing is off.
+
+procedure good_pixels (im, imfit, gl, sf)
+
+pointer	im	# input image
+pointer	imfit	# pointer to imsurfit header structure
+pointer	gl	# pointer to good region list
+pointer	sf	# pointer to the surface descriptor
+
+int	lp, lineno, prevlineno, ncols, nlines, npts, nranges, ier, ijunk
+int	max_nranges
+pointer	sp, colsfit, lines, buf, fbuf, wgt, ranges
+int	prl_nextlineno(), prl_eqlines(), prl_get_ranges(), is_expand_ranges()
+int	is_choose_rangesr()
+pointer	imgl2r()
+
+errchk	smark, salloc, sfree, imgl2r
+errchk	isinit, islfit, islrefit, issolve
+errchk	prl_nextlineno, prl_eqlines, prl_get_ranges
+errchk	is_choose_rangesr
+
+begin
+	ncols = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	max_nranges = ncols
+
+	# Initialize the surface fit.
+	call isinit (sf, SURFACE_TYPE(imfit), XORDER(imfit), YORDER(imfit),
+	    CROSS_TERMS(imfit), ncols, nlines)
+
+	# Allocate temporary space for fitting.
+	call smark (sp)
+	call salloc (colsfit, ncols, TY_INT)
+	call salloc (lines, nlines, TY_INT)
+	call salloc (fbuf, ncols, TY_REAL)
+	call salloc (wgt, ncols, TY_REAL)
+	call salloc (ranges, 3 * max_nranges + 1, TY_INT)
+	call amovkr (1., Memr[wgt], ncols)
+
+	# Intialize counters and pointers.
+	lp = 0
+	lineno = 0
+	prevlineno = 0
+
+	# Loop over those lines to be fit.
+	while (prl_nextlineno (gl, lineno) != EOF) {
+
+	    # Read in the image line.
+	    buf = imgl2r (im, lineno)
+	    if (buf == EOF)
+		call error (0, "GOOD_PIXELS: Error reading image.")
+
+	    # Get the ranges for that image line.
+	    nranges = prl_get_ranges (gl, lineno, Memi[ranges], max_nranges)
+	    if (nranges == 0)
+		next
+
+	    # If ranges are not equal to previous line fit else refit.
+	    if (lp == 0 || prl_eqlines (gl, lineno, prevlineno) == NO) {
+		npts = is_expand_ranges (Memi[ranges], Memi[colsfit], ncols)
+		ijunk = is_choose_rangesr (Memi[colsfit], Memr[buf], Memr[fbuf],
+		    npts, 1, ncols)
+		call islfit (sf, Memi[colsfit], lineno, Memr[fbuf], Memr[wgt],
+		    npts, SF_USER, ier)
+	    } else {
+		ijunk = is_choose_rangesr (Memi[colsfit], Memr[buf], Memr[fbuf],
+		    npts, 1, ncols)
+		call islrefit (sf, Memi[colsfit], lineno, Memr[fbuf], Memr[wgt])
+	    }
+
+	    # Handle fitting errors.
+	    switch (ier) {
+	    case NO_DEG_FREEDOM:
+		call eprintf ("Warning: Too few columns to fit line: %d\n")
+		    call pargi (lineno)
+	    case SINGULAR:
+		call eprintf ("Warning: Solution singular for line: %d\n")
+		    call pargi (lineno)
+		Memi[lines+lp] = lineno
+		lp = lp + 1
+	    default:
+		Memi[lines+lp] = lineno
+		lp = lp + 1
+	    }
+
+	    prevlineno = lineno
+	}
+
+	# Solve the surface.
+	call issolve (sf, Memi[lines], lp, ier)
+
+	# Handle fitting errors.
+	switch (ier) {
+	case NO_DEG_FREEDOM:
+	    call error (0, "GOOD_PIXELS: Cannot fit surface.")
+	case SINGULAR:
+	    call eprintf ("Warning: Solution singular for surface.\n")
+	default:
+	    # everything OK
+	}
+
+	# Free space.
+	call sfree (sp)
+end
+
+
+# ALL_MEDIANS -- Get surface when median processor on and all
+# pixels good.
+
+procedure  all_medians (im, imfit, sf)
+
+pointer	im	# input image
+pointer	imfit	# pointer to the imsurfit header structure
+pointer	sf	# pointer to the surface descriptor
+
+int	i, lp, cp, op, lineno, x1, x2, y1, y2, ier
+int	nimcols, nimlines, ncols, nlines, npts
+pointer	sp, cols, lines, wgt, z, med, sbuf, lbuf, buf
+
+pointer	imgs2r()
+real	asokr()
+errchk	salloc, sfree, smark
+errchk	isinit, islfit, islrefit, issolve
+
+begin
+	# Determine the number of lines and columns for a median processed
+	# image.
+	nimcols = IM_LEN(im,1)
+	if (mod (int (IM_LEN(im,1)), XMEDIAN(imfit)) != 0)
+	    ncols = IM_LEN(im,1) / XMEDIAN(imfit) + 1
+	else
+	    ncols = IM_LEN(im,1) / XMEDIAN(imfit)
+	nimlines = IM_LEN(im,2)
+	if (mod (int (IM_LEN(im,2)), YMEDIAN(imfit)) != 0)
+	    nlines = IM_LEN(im,2) / YMEDIAN(imfit) + 1
+	else
+	    nlines = IM_LEN(im,2) / YMEDIAN(imfit)
+
+	# Initialize the surface fitting.
+	call isinit (sf, SURFACE_TYPE(imfit), XORDER(imfit), YORDER(imfit),
+	   CROSS_TERMS(imfit), ncols, nlines)
+
+	# Allocate workin memory.
+	call smark (sp)
+	call salloc (cols, ncols, TY_INT)
+	call salloc (wgt, ncols, TY_REAL)
+	call salloc (lines, nlines, TY_INT)
+	call salloc (z, ncols, TY_REAL)
+	call salloc (med, XMEDIAN(imfit) * YMEDIAN(imfit), TY_REAL)
+
+	# Intialize the x and weight arrays.
+	do i = 1, ncols
+	    Memi[cols - 1 + i] = i
+	call amovkr (1.0, Memr[wgt], ncols)
+
+	# Loop over image sections.
+	lp = 0
+	lineno = 1
+	for (y1 = 1; y1 <= nimlines; y1 = y1 + YMEDIAN(imfit)) {
+
+	    # Get image section.
+	    y2 = min (y1 + YMEDIAN(imfit) - 1, nimlines)
+	    sbuf = imgs2r (im, 1, nimcols, y1, y2)
+	    if (sbuf == EOF)
+		call error (0, "Error reading image section.")
+
+	    # Loop over median boxes.
+	    cp = 0
+	    for (x1 = 1; x1 <= nimcols; x1 = x1 + XMEDIAN(imfit)) {
+
+		x2 = min (x1 + XMEDIAN(imfit) - 1, nimcols)
+		npts = x2 - x1 + 1
+		lbuf = sbuf - 1 + x1
+
+		# Loop over lines in the median box.
+		op = 0
+		buf = lbuf
+		for (i = 1; i <= y2 - y1 + 1; i = i + 1) {
+		   call amovr (Memr[buf], Memr[med+op], npts)
+		   op = op + npts
+		   buf = buf + nimcols
+		}
+
+		# Calculate the median.
+		Memr[z+cp] = asokr (Memr[med], op, (op + 1) / 2)
+		cp = cp + 1
+
+	    }
+
+	    # Fit each image "line".
+	    if (y1 == 1)
+		call islfit (sf, Memi[cols], lineno, Memr[z], Memr[wgt],
+		     ncols, SF_USER, ier)
+	    else
+		call islrefit (sf, Memi[cols], lineno, Memr[z], Memr[wgt])
+
+	    # Handle fitting errors.
+	    switch (ier) {
+	    case NO_DEG_FREEDOM:
+		call eprintf ("Warning: Too few columns to fit line: %d\n")
+		    call pargi (lineno)
+	    case SINGULAR:
+		call eprintf ("Warning: Solution singular for line: %d\n")
+		    call pargi (lineno)
+		Memi[lines + lp] = lineno
+		lp = lp + 1
+	    default:
+		Memi[lines + lp] = lineno
+		lp = lp + 1
+	    }
+
+	    lineno = lineno + 1
+	}
+
+	# Solve th surface.
+	call issolve (sf, Memi[lines], lp, ier)
+
+	# Handle fitting errors.
+	switch (ier) {
+	case NO_DEG_FREEDOM:
+	    call error (0, "ALL_MEDIANS: Cannot fit surface.")
+	case SINGULAR:
+	    call eprintf ("Warning: Solution singular for surface.\n")
+	default:
+	    # everything OK
+	}
+
+	# Free space
+	call sfree (sp)
+end
+
+
+# GOOD_MEDIANS -- Procedure to fetch medians when the good regions
+# list is defined.
+
+procedure good_medians (im, imfit, gl, sf)
+
+pointer	im	# input image
+pointer	imfit	# pointer to surface descriptor structure
+pointer	gl	# pointer to good regions list
+pointer	sf	# pointer the surface descriptor
+
+int	i, cp, lp, x1, x2, y1, y2, ier, ntemp
+int	nimcols, nimlines, ncols, nlines, nranges, nbox, nxpts
+int	lineno, current_line, lines_per_box, max_nranges
+pointer	sp, colsfit, cols, lines, wgt, npts, lbuf, med, mbuf, z, ranges
+
+int	prl_get_ranges(), prl_nextlineno(), is_expand_ranges()
+int	is_choose_rangesr()
+pointer	imgl2r()
+real	asokr()
+errchk	smark, salloc, sfree, imgl2r
+errchk	isinit, islfit, issolve
+errchk	prl_get_ranges, prl_nextlineno, is_choose_rangesr()
+
+begin
+	# Determine the number of lines and columns for a median processed
+	# image.
+	nimcols = IM_LEN(im,1)
+	if (mod (int (IM_LEN(im,1)), XMEDIAN(imfit)) != 0)
+	    ncols = IM_LEN(im,1) / XMEDIAN(imfit) + 1
+	else
+	    ncols = IM_LEN(im,1) / XMEDIAN(imfit)
+	nimlines = IM_LEN(im,2)
+	if (mod (int (IM_LEN(im,2)), YMEDIAN(imfit)) != 0)
+	    nlines = IM_LEN(im,2) / YMEDIAN(imfit) + 1
+	else
+	    nlines = IM_LEN(im,2) / YMEDIAN(imfit)
+	nbox = XMEDIAN(imfit) * YMEDIAN(imfit)
+	max_nranges = nimcols
+
+	# Initialize the surface fitting.
+	call isinit (sf, SURFACE_TYPE(imfit), XORDER(imfit), YORDER(imfit),
+	    CROSS_TERMS(imfit), ncols, nlines)
+
+	# Allocate working memory.
+	call smark (sp)
+	call salloc (colsfit, nimcols, TY_INT)
+	call salloc (cols, ncols, TY_INT)
+	call salloc (npts, ncols, TY_INT)
+	call salloc (lines, nlines, TY_INT)
+	call salloc (wgt, ncols, TY_REAL)
+	call salloc (med, nbox * ncols, TY_REAL)
+	call salloc (z, ncols, TY_REAL)
+	call salloc (ranges, 3 * max_nranges + 1, TY_INT)
+	call amovkr (1., Memr[wgt], ncols)
+
+	# Loop over median boxes in y.
+	lp = 0
+	lineno = 0
+	for (y1 = 1; y1 <= nimlines; y1 = y1 + YMEDIAN(imfit)) {
+
+	    lineno = lineno + 1
+	    current_line = y1 - 1
+	    y2 = min (y1 + YMEDIAN(imfit) - 1, nimlines)
+
+	    # If lines not in range, next image section.
+	    lines_per_box = 0
+	    while (prl_nextlineno (gl, current_line) != EOF) {
+		if (current_line > y2)
+		    break
+		lines_per_box = lines_per_box + 1
+	    }
+	    if (lines_per_box < (YMEDIAN(imfit) * (MEDIAN_PERCENT(imfit)/100.)))
+		next
+
+	    # Loop over the image lines.
+	    call aclri (Memi[npts], ncols)
+	    do i = y1, y2 {
+
+		# Get image line, and check the good regions list.
+		lbuf = imgl2r (im, i)
+		nranges = prl_get_ranges (gl, i, Memi[ranges], max_nranges)
+		if (nranges == 0)
+		    next
+		nxpts = is_expand_ranges (Memi[ranges], Memi[colsfit], nimcols)
+
+		# Loop over the median boxes in x.
+		cp= 0
+		mbuf = med
+		for (x1 = 1; x1 <= nimcols; x1 = x1 + XMEDIAN(imfit)) {
+		    x2 = min (x1 + XMEDIAN(imfit) - 1, nimcols)
+		    ntemp = is_choose_rangesr (Memi[colsfit], Memr[lbuf],
+		        Memr[mbuf+Memi[npts+cp]], nxpts, x1, x2)
+		    Memi[npts+cp] = Memi[npts+cp] + ntemp
+		    mbuf = mbuf + nbox
+		    cp = cp + 1
+		}
+	    }
+
+	    # Calculate the medians.
+	    nxpts = 0
+	    mbuf = med
+	    do i = 1, ncols {
+		if (Memi[npts+i-1] > ((MEDIAN_PERCENT(imfit) / 100.) * nbox)) {
+		    Memr[z+nxpts] = asokr (Memr[mbuf], Memi[npts+i-1],
+			(Memi[npts+i-1] + 1) / 2)
+		    Memi[cols+nxpts] = i
+		    nxpts = nxpts + 1
+		}
+		mbuf = mbuf + nbox
+	    }
+
+	    # Fit the line.
+	    call islfit (sf, Memi[cols], lineno, Memr[z], Memr[wgt], nxpts,
+		SF_USER, ier)
+
+	    # Handle fitting errors.
+	    switch (ier) {
+	    case NO_DEG_FREEDOM:
+		call eprintf ("Warning: Too few columns to fit line: %d\n")
+		    call pargi (lineno)
+	    case SINGULAR:
+		call eprintf ("Warning: Solution singular for line: %d\n")
+		    call pargi (lineno)
+		Memi[lines+lp] = lineno
+		lp = lp + 1
+	    default:
+		Memi[lines+lp] = lineno
+		lp = lp + 1
+	    }
+	}
+
+	# Solve the surface.
+	call issolve (sf, Memi[lines], lp, ier)
+
+	# Handle fitting errors.
+	switch (ier) {
+	case NO_DEG_FREEDOM:
+	    call error (0, "GOOD_MEDIANS: Cannot fit surface.")
+	case SINGULAR:
+	    call eprintf ("Warning: Solution singular for surface.")
+	default:
+	    # everyting OK
+	}
+
+	# Free space.
+	call sfree (sp)
+end
+
+
+# SET_OUTIMAGE -- Procedure to write an output image of the desired type.
+
+procedure set_outimage (imin, imout, imfit, sf, rl)
+
+pointer	imin	# input image
+pointer	imout	# output image
+pointer	imfit	# pointer to the imsurfut header structure
+pointer	sf	# pointer to the surface descriptor
+pointer	rl	# pointer to the rejected pixel list regions list
+
+int	i, k, ncols, nlines, max_nranges
+long	u[IM_MAXDIM], v[IM_MAXDIM]
+real	b1x, b2x, b1y, b2y
+pointer	sp, x, y, inbuf, outbuf, ranges
+
+int	impnlr(), imgnlr()
+real	ims_divzero()
+extern	ims_divzero
+errchk	malloc, mfree, imgnlr, impnlr
+
+begin
+	ncols = IM_LEN(imin,1)
+	nlines = IM_LEN(imin,2)
+	max_nranges = ncols
+
+	# Calculate transformation constants from real coordinates to
+	# median coordinates if median processing specified. 
+	if (MEDIAN(imfit) == YES) {
+	    b1x = (1. + XMEDIAN(imfit)) / (2. * XMEDIAN(imfit))
+	    b2x = (2. * ncols + XMEDIAN(imfit) - 1.) / (2. * XMEDIAN(imfit))
+	    b1y = (1. + YMEDIAN(imfit)) / (2. * YMEDIAN(imfit))
+	    b2y = (2. * nlines + YMEDIAN(imfit) - 1.) / (2. * YMEDIAN(imfit))
+	}
+
+	# Allocate space for x coordinates, initialize to image coordinates
+	# and transform to median coordinates.
+	call smark (sp)
+	call salloc (x, ncols, TY_REAL)
+	call salloc (y, ncols, TY_REAL)
+	call salloc (ranges, 3 * max_nranges + 1, TY_INT)
+
+	# Intialize the x array.
+	do i = 1, ncols
+	    Memr[x - 1 + i] = i
+	if (MEDIAN(imfit) == YES)
+	    call amapr (Memr[x], Memr[x], ncols, 1., real (ncols), b1x, b2x)
+
+	# loop over the images lines
+	call amovkl (long (1), v, IM_MAXDIM)
+	call amovkl (long (1), u, IM_MAXDIM)
+	do i = 1, nlines {
+
+	    # Get input and output image buffers.
+	    if (TYPE_OUTPUT(imfit) != FIT) {
+		if (imgnlr (imin, inbuf, v) == EOF)
+		    call error (0, "Error reading input image.")
+	    }
+	    if (impnlr (imout, outbuf, u) == EOF)
+		call error (0, "Error writing output image.")
+
+	    # Intialize y coordinates to image coordinates, and
+	    # transform to median coordinates.
+	    if (MEDIAN(imfit) == YES) {
+	        Memr[y] = real (i)
+		call amapr (Memr[y], Memr[y], 1, 1., real (nlines),
+			    b1y, b2y)
+		call amovkr (Memr[y], Memr[y+1], (ncols-1))
+	    } else
+		call amovkr (real (i), Memr[y], ncols)
+
+	    # Write output image.
+	    switch (TYPE_OUTPUT(imfit)) {
+	    case FIT:
+	        call isvector (sf, Memr[x], Memr[y], Memr[outbuf], ncols)
+	    case CLEAN:
+		call clean_line (Memr[x], Memr[y], Memr[inbuf], ncols, nlines,
+		   rl, sf, i, NGROW(imfit)) 
+		call amovr (Memr[inbuf], Memr[outbuf], ncols)
+	    case RESID:
+	        call isvector (sf, Memr[x], Memr[y], Memr[outbuf], ncols)
+		call asubr (Memr[inbuf], Memr[outbuf], Memr[outbuf], ncols)
+	    case RESP:
+	        call isvector (sf, Memr[x], Memr[y], Memr[outbuf], ncols)
+		if (IS_INDEF(DIV_MIN(imfit))) {
+		    iferr (call adivr (Memr[inbuf], Memr[outbuf], Memr[outbuf],
+		        ncols))
+			call advzr (Memr[inbuf], Memr[outbuf], Memr[outbuf],
+			    ncols, ims_divzero)
+		} else {
+		    do k = 1, ncols {
+			if (Memr[outbuf-1+k] < DIV_MIN(imfit))
+			    Memr[outbuf-1+k] = 1.
+			else
+			    Memr[outbuf-1+k] = Memr[inbuf-1+k] /
+			    Memr[outbuf-1+k]
+		    }
+		}
+	    default:
+		call error (0, "SET_OUTIMAGE: Unknown output type.")
+	    }
+	}
+
+	# Free space
+	call sfree (sp)
+end
+
+
+# CLEAN_LINE -- Procedure to set weights of rejected points to zero
+
+procedure clean_line (x, y, z, ncols, nlines, rl, sf, line, ngrow)
+
+real	x[ARB]		# array of weights set to 1
+real	y		# y value of line
+real	z[ARB]		# line of data
+int	ncols		# number of image columns
+int	nlines		# number of image lines
+pointer	rl		# pointer to reject pixel list
+pointer	sf		# surface fitting
+int	line		# line number
+int	ngrow		# radius for region growing
+
+int	cp, j, k, nranges, dist, yreg_min, yreg_max, xreg_min, xreg_max
+pointer	sp, branges
+real	r2
+int	prl_get_ranges(), is_next_number()
+real	iseval()
+
+begin
+	call smark (sp)
+	call salloc (branges, 3 * ncols + 1, TY_INT)
+
+	r2 = ngrow ** 2
+	yreg_min = max (1, line - ngrow)
+	yreg_max = min (nlines, line + ngrow)
+
+	do j = yreg_min, yreg_max {
+	    nranges = prl_get_ranges (rl, j, Memi[branges], ncols)
+	    if (nranges == 0)
+		next
+	    dist = int (sqrt (r2 - (j - line) ** 2))
+	    cp = 0
+	    while (is_next_number (Memi[branges], cp) != EOF) {
+		xreg_min = max (1, cp - dist)
+		xreg_max = min (ncols, cp + dist)
+		do k = xreg_min, xreg_max
+		    z[k] = iseval (sf, x[k], y)
+		cp = xreg_max
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# DO_REJECT -- Procedure to detect rejected pixels in an image.
+
+procedure do_reject (im, imfit, gl, sf, rl)
+
+pointer		im	# pointer to in put image
+pointer		imfit	# pointer to image fitting structure
+pointer		gl	# pointer to good regions list
+pointer		sf	# pointer to surface descriptor
+pointer		rl	# pointer to rejected pixel list
+
+int	niter, nrejects
+real	sigma
+int	detect_rejects()
+real	get_sigma()
+errchk  prl_init, detect_rejects, get_sigma, refit_surface
+
+begin
+	# Initialize rejected pixel list.
+	call prl_init (rl, int(IM_LEN(im,1)), int(IM_LEN(im,2)))
+
+	# Do an iterative rejection cycle on the image.
+	niter = 0
+	repeat {
+
+	    # Get the sigma of the fit.
+	    sigma = get_sigma (im, gl, sf, rl)
+
+	    # Detect rejected pixels.
+	    nrejects = detect_rejects (im, imfit, gl, sf, rl, sigma)
+
+	    # If no rejected pixels quit, else refit surface.
+	    if (nrejects == 0 || NITER(imfit) == 0)
+		break
+	    call refit_surface (im, imfit, gl, sf, rl)
+
+	    niter = niter + 1
+
+	} until (niter == NITER(imfit))
+end
+
+
+# REFIT_SURFACE -- Procedure tp refit the surface.
+
+procedure refit_surface (im, imfit, gl, sf, rl)
+
+pointer		im	# pointer to image
+pointer		imfit	# pointer to surface fitting structure
+pointer		gl	# pointer to good regions list
+pointer		sf	# pointer to surface descriptor
+pointer		rl	# pointer to rejected pixels list
+
+int	i, ijunk, lp, ier, max_nranges
+int	ncols, nlines, npts, nfree, nrejects, nranges, ncoeff
+pointer	sp, cols, colsfit, lines, buf, fbuf, wgt, granges
+
+int	prl_get_ranges(), grow_regions(), is_expand_ranges()
+int	is_choose_rangesr()
+pointer	imgl2r()
+errchk	smark, salloc, sfree, imgl2r
+errchk	iscoeff, islfit, issolve
+errchk	prl_get_ranges, grow_regions
+errchk	is_choose_rangesr
+
+begin
+	ncols = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	max_nranges = ncols
+
+	# Allocate up temporary storage.
+	call smark (sp)
+	call salloc (cols, ncols, TY_INT)
+	call salloc (colsfit, ncols, TY_INT)
+	call salloc (lines, nlines, TY_INT)
+	call salloc (fbuf, ncols, TY_REAL)
+	call salloc (wgt, ncols, TY_REAL)
+	call salloc (granges, 3 * max_nranges + 1, TY_INT)
+
+	# Initialize columns.
+	do i = 1, ncols
+	    Memi[cols+i-1] = i
+	call amovi (Memi[cols], Memi[colsfit], ncols)
+
+	# Get number of coefficients.
+	switch (SURFACE_TYPE(imfit)) {
+	case SF_LEGENDRE, SF_CHEBYSHEV:
+	    ncoeff = XORDER(imfit)
+	case SF_SPLINE3:
+	    ncoeff = XORDER(imfit) + 3
+	case SF_SPLINE1:
+	    ncoeff = XORDER(imfit) + 1
+	}
+
+	# Refit affected lines and solve for surface.
+	lp = 0
+	do i = 1, nlines {
+
+	    # Determine whether image line is good.
+	    if (gl != NULL) {
+		nranges = prl_get_ranges (gl, i, Memi[granges], max_nranges)
+		if (nranges == 0)
+		    next
+	    }
+
+	    # Define rejected points with region growing.
+	    call amovkr (1., Memr[wgt], ncols)
+	    nrejects = grow_regions (Memr[wgt], ncols, nlines, rl, i,
+	        NGROW(imfit))
+
+	    # Get number of data points.
+	    if (gl == NULL)
+		npts = ncols
+	    else
+		npts = is_expand_ranges (Memi[granges], Memi[colsfit], ncols)
+	    nfree = npts - nrejects
+
+	    # If no rejected pixels skip to next line.
+	    if (nrejects == 0) {
+		if (nfree >= ncoeff ) {
+		    Memi[lines+lp] = i
+		    lp = lp + 1
+		}
+		next
+	    }
+
+	    # Read in image line.
+	    buf = imgl2r (im, i)
+	    if (buf == EOF)
+		call error (0, "REFIT_SURFACE: Error reading image.")
+
+	    # Select the data.
+	    if (gl == NULL) {
+		npts = ncols
+	        if (nfree >= ncoeff)
+	            call islfit (sf, Memi[colsfit], i, Memr[buf], Memr[wgt],
+		        npts, SF_USER, ier)
+	        else
+		    ier = NO_DEG_FREEDOM
+	    } else {
+		ijunk = is_choose_rangesr (Memi[colsfit], Memr[buf],
+		    Memr[fbuf], npts, 1, ncols)
+		ijunk = is_choose_rangesr (Memi[colsfit], Memr[wgt], Memr[wgt],
+		    npts, 1, ncols)
+	        if (nfree >= ncoeff)
+	            call islfit (sf, Memi[colsfit], i, Memr[fbuf], Memr[wgt],
+		        npts, SF_USER, ier)
+	        else
+		    ier = NO_DEG_FREEDOM
+	    }
+
+	    # Evaluate fitting errors.
+	    switch (ier) {
+	    case NO_DEG_FREEDOM:
+		call eprintf ("REFIT_SURFACE: Too few points to fit line: %d\n")
+		    call pargi (i)
+	    case SINGULAR:
+		call eprintf ("REFIT_SURFACE: Solution singular for line: %d\n")
+		    call pargi (i)
+		Memi[lines+lp] = i
+		lp = lp + 1
+	    default:
+		Memi[lines+lp] = i
+		lp = lp + 1
+	    }
+	}
+
+	# Resolve surface.
+	call issolve (sf, Memi[lines], lp, ier)
+
+	# Evaluate fitting errors for surface.
+	switch (ier) {
+	case NO_DEG_FREEDOM:
+	    call error (0, "REFIT_SURFACE: Too few points to fit surface\n")
+	case SINGULAR:
+	    call eprintf ("REFIT_SURFACE: Solution singular for surface\n")
+	default:
+	    # everything OK
+	}
+
+	call sfree (sp)
+
+end
+
+
+# GROW_REGIONS -- Procedure to set weights of rejected points to zero.
+
+int procedure grow_regions (wgt, ncols, nlines, rl, line, ngrow)
+
+real	wgt[ARB]	# array of weights set to 1
+int	ncols		# number of image columnspoints
+int	nlines		# number of images lines
+pointer	rl		# pointer to reject pixel list
+int	line		# line number
+int	ngrow		# radius for region growing
+
+int	cp, j, k, nrejects, nranges, max_nranges
+int	dist, yreg_min, yreg_max, xreg_min, xreg_max
+pointer	sp, branges
+real	r2
+int	prl_get_ranges(), is_next_number()
+errchk	smark, salloc, sfree
+errchk	prl_get_ranges, is_next_number 
+
+begin
+	max_nranges = ncols
+
+	call smark (sp)
+	call salloc (branges, 3 * max_nranges + 1, TY_INT)
+
+	r2 = ngrow ** 2
+	nrejects = 0
+	yreg_min = max (1, line - ngrow)
+	yreg_max = min (nlines, line + ngrow)
+
+	do j = yreg_min, yreg_max {
+	    nranges = prl_get_ranges (rl, j, Memi[branges], max_nranges)
+	    if (nranges == 0)
+		next
+	    dist = int (sqrt (r2 - (j - line) ** 2))
+	    cp = 0
+	    while (is_next_number (Memi[branges], cp) != EOF) {
+		xreg_min = max (1, cp - dist)
+		xreg_max = min (ncols, cp + dist)
+		do k = xreg_min, xreg_max {
+		    if (wgt[k] > 0.) {
+			wgt[k] = 0.
+			nrejects = nrejects + 1
+		    }
+		}
+		cp = xreg_max
+	    }
+	}
+
+	call sfree (sp)
+	return (nrejects)
+end
+
+
+# GET_SIGMA -- Procedure to calculate the sigma of the surface fit
+
+real procedure get_sigma (im, gl, sf, rl)
+
+pointer		im	# pointer to image
+pointer		gl	# pointer to good pixel list
+pointer		sf	# pointer to surface deascriptor
+pointer		rl	# pointer to rejected pixel list
+
+int	i, ijunk, cp, nranges, npts, ntpts, ncols, nlines, max_nranges
+pointer	sp, colsfit, x, xfit, y, zfit, buf, fbuf, wgt, granges, branges
+real	sum, sigma
+int	prl_get_ranges(), is_next_number(), is_expand_ranges()
+int	is_choose_rangesr()
+pointer	imgl2r()
+real	asumr(), awssqr()
+errchk	smark, salloc, sfree, imgl2r
+
+begin
+	ncols = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	max_nranges = ncols
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (colsfit, ncols, TY_INT)
+	call salloc (x, ncols, TY_REAL)
+	call salloc (xfit, ncols, TY_REAL)
+	call salloc (y, ncols, TY_REAL)
+	call salloc (fbuf, ncols, TY_REAL)
+	call salloc (zfit, ncols, TY_REAL)
+	call salloc (wgt, ncols, TY_REAL)
+	call salloc (granges, 3 * max_nranges + 1, TY_INT)
+	call salloc (branges, 3 * max_nranges + 1, TY_INT)
+
+	# Intialize the x array.
+	do i = 1, ncols
+	    Memr[x+i-1] = i
+	call amovr (Memr[x], Memr[xfit], ncols)
+
+	sum = 0.
+	sigma = 0.
+	ntpts = 0
+
+	# Loop over the image.
+	do i = 1, nlines {
+
+	    # Check that line is in range.
+	    if (gl != NULL) {
+		nranges = prl_get_ranges (gl, i, Memi[granges], max_nranges)
+		if (nranges == 0)
+		    next
+		npts = is_expand_ranges (Memi[granges], Memi[colsfit], ncols)
+	    }
+	    
+	    # Read in image.
+	    buf = imgl2r (im, i)
+	    if (buf == EOF)
+		call error (0, "GET_SIGMA: Error reading image.")
+
+	    # Select appropriate data and fit.
+	    call amovkr (real (i), Memr[y], ncols)
+	    if (gl == NULL) {
+		npts = ncols
+	        call isvector (sf, Memr[xfit], Memr[y], Memr[zfit], npts)
+	        call asubr (Memr[buf], Memr[zfit], Memr[zfit], npts)
+	    } else {
+		ijunk = is_choose_rangesr (Memi[colsfit], Memr[x], Memr[xfit],
+		    npts, 1, ncols)
+		ijunk = is_choose_rangesr (Memi[colsfit], Memr[buf], Memr[fbuf],
+		    npts, 1, ncols)
+	        call isvector (sf, Memr[xfit], Memr[y], Memr[zfit], npts)
+	        call asubr (Memr[fbuf], Memr[zfit], Memr[zfit], npts)
+	    }
+
+	    # Get ranges of rejected pixels for the line and set weights.
+	    call amovkr (1., Memr[wgt], ncols)
+	    nranges = prl_get_ranges (rl, i, Memi[branges], max_nranges)
+	    if (nranges > 0) {
+		cp = 0
+		while (is_next_number (Memi[branges], cp) != EOF)
+		    Memr[wgt+cp-1] = 0.
+		if (gl != NULL)
+		    ijunk = is_choose_rangesr (Memi[colsfit], Memr[wgt],
+		        Memr[wgt], npts, 1, ncols)
+	    }
+
+	    # Calculate sigma.
+	    sigma = sigma + awssqr (Memr[zfit], Memr[wgt], npts)
+	    ntpts = ntpts + asumr (Memr[wgt], npts)
+	}
+
+	call sfree (sp)
+
+	return (sqrt (sigma / (ntpts - 1)))
+end
+
+
+# DETECT_REJECTS - Procedure to detect rejected pixels.
+
+int procedure detect_rejects (im, imfit, gl, sf, rl, sigma)
+
+pointer		im	# pointer to image
+pointer		imfit	# pointer to surface fitting structure
+pointer		gl	# pointer to good pixel list
+pointer		sf	# pointer to surface descriptor
+pointer		rl	# pointer to rejected pixel list
+real		sigma	# standard deviation of fit
+
+int	i, j, ijunk, cp, ncols, nlines, npts, nranges, nlrejects, ntrejects
+int	norejects, max_nranges
+pointer	sp, granges, branges, x, xfit, cols, colsfit, y, zfit, buf, fbuf
+pointer	wgt, list
+real	upper, lower
+
+int	prl_get_ranges(), is_next_number(), is_make_ranges(), is_expand_ranges()
+int	is_choose_rangesr()
+pointer	imgl2r()
+
+begin
+	ncols = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	max_nranges = ncols
+
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (x, ncols, TY_REAL)
+	call salloc (xfit, ncols, TY_REAL)
+	call salloc (cols, ncols, TY_INT)
+	call salloc (colsfit, ncols, TY_INT)
+	call salloc (y, ncols, TY_REAL)
+	call salloc (fbuf, ncols, TY_REAL)
+	call salloc (zfit, ncols, TY_REAL)
+	call salloc (wgt, ncols, TY_REAL)
+	call salloc (granges, 3 * max_nranges + 1, TY_INT)
+	call salloc (branges, 3 * max_nranges + 1, TY_INT)
+	call salloc (list, ncols, TY_INT)
+
+	# Intialize x and column values.
+	do i = 1, ncols {
+	    Memi[cols+i-1] = i
+	    Memr[x+i-1] = i
+	}
+	call amovr (Memr[x], Memr[xfit], ncols)
+	call amovi (Memi[cols], Memi[colsfit], ncols)
+
+	ntrejects = 0
+	if (LOWER(imfit) <= 0.0)
+	    lower = -MAX_REAL
+	else
+	    lower = -sigma * LOWER(imfit)
+	if (UPPER(imfit) <= 0.0)
+	    upper = MAX_REAL
+	else
+	    upper = sigma * UPPER(imfit)
+
+	# Loop over the image.
+	do i = 1, nlines {
+
+	    # Get ranges if appropriate.
+	    if (gl != NULL) {
+		nranges = prl_get_ranges (gl, i, Memi[granges], max_nranges)
+		if (nranges == 0)
+		    next
+		npts = is_expand_ranges (Memi[granges], Memi[colsfit], ncols) 
+	    }
+
+	    # Read in image.
+	    buf = imgl2r (im, i)
+	    if (buf == EOF)
+		call error (0, "GET_SIGMA: Error reading image.")
+
+	    # Select appropriate data and fit.
+	    call amovkr (real (i), Memr[y], ncols)
+	    if (gl == NULL) {
+		npts = ncols
+	        call isvector (sf, Memr[xfit], Memr[y], Memr[zfit], npts)
+	        call asubr (Memr[buf], Memr[zfit], Memr[zfit], npts)
+	    } else {
+		ijunk = is_choose_rangesr (Memi[colsfit], Memr[x], Memr[xfit],
+		    npts, 1, ncols)
+		ijunk = is_choose_rangesr (Memi[colsfit], Memr[buf], Memr[fbuf],
+		    npts, 1, ncols)
+	        call isvector (sf, Memr[xfit], Memr[y], Memr[zfit], npts)
+	        call asubr (Memr[fbuf], Memr[zfit], Memr[zfit], npts)
+	    }
+
+	    # Get ranges of rejected pixels for the line and set weights.
+	    call amovkr (1., Memr[wgt], ncols)
+	    nranges = prl_get_ranges (rl, i, Memi[branges], max_nranges)
+	    norejects = 0
+	    if (nranges > 0) {
+		cp = 0
+		while (is_next_number (Memi[branges], cp) != EOF) {
+		    Memi[list+norejects] = cp
+		    norejects = norejects + 1
+		    Memr[wgt+cp-1] = 0.
+		}
+		if (gl != NULL)
+		    ijunk = is_choose_rangesr (Memi[colsfit], Memr[wgt],
+		        Memr[wgt], npts, 1, ncols)
+	    }
+
+	    # Detect deviant pixels.
+	    nlrejects = 0
+	    do j = 1, npts {
+		if ((Memr[zfit+j-1] < lower || Memr[zfit+j-1] > upper) &&
+		    Memr[wgt+j-1] != 0.) {
+		    Memi[list+norejects+nlrejects] = Memi[colsfit+j-1]
+		    nlrejects = nlrejects + 1
+		}
+	    }
+
+	    # Add to rejected pixel list.
+	    if (nlrejects > 0) {
+		call asrti (Memi[list], Memi[list], norejects + nlrejects)
+		nranges = is_make_ranges (Memi[list], norejects + nlrejects,
+		    Memi[granges], max_nranges)
+		call prl_put_ranges (rl, i, i, Memi[granges])
+	    }
+
+	    ntrejects = ntrejects + nlrejects
+	}
+
+	call sfree (sp)
+	return (ntrejects)
+end
+
+
+# AWSSQR -- Procedure to calculate the weighted sum of the squares
+
+real procedure awssqr (a, w, npts)
+
+real	a[npts]		# array of data
+real	w[npts]		# array of points
+int	npts		# number of data points
+
+int	i
+real	sum
+
+begin
+	sum = 0.
+	do i = 1, npts
+	    sum = sum + w[i] * a[i] ** 2
+
+	return (sum)
+end
+
+
+# IMS_DIVZER0 -- Return 1. on a divide by zero
+
+real	procedure ims_divzero (x)
+
+real	x
+
+begin
+	return (1.)
+end
diff --git a/pkg/images/imfit/src/mkpkg b/pkg/images/imfit/src/mkpkg
new file mode 100644
index 00000000..3bc27b8f
--- /dev/null
+++ b/pkg/images/imfit/src/mkpkg
@@ -0,0 +1,15 @@
+# Library for the IMAGES IMFIT Subpackage Tasks
+
+$checkout libpkg.a ../../
+$update   libpkg.a
+$checkin  libpkg.a ../../
+$exit
+
+libpkg.a:
+	fit1d.x		<imhdr.h> <pkg/gtools.h> <error.h>
+	imsurfit.x	imsurfit.h <math/surfit.h> <imhdr.h> <mach.h>
+	pixlist.x       pixlist.h
+	ranges.x        <ctype.h> <mach.h>
+	t_imsurfit.x	imsurfit.h <error.h> <imhdr.h>
+	t_lineclean.x	<imhdr.h> <pkg/gtools.h>
+	;
diff --git a/pkg/images/imfit/src/pixlist.h b/pkg/images/imfit/src/pixlist.h
new file mode 100644
index 00000000..39c875a9
--- /dev/null
+++ b/pkg/images/imfit/src/pixlist.h
@@ -0,0 +1,11 @@
+# PIXEL LIST descriptor structure
+
+define	LEN_PLSTRUCT	10
+
+define	PRL_NCOLS		Memi[$1]	# number of columns
+define	PRL_NLINES		Memi[$1+1]	# number of lines
+define	PRL_LINES		Memi[$1+2]	# pointer to the line offsets
+define	PRL_LIST		Memi[$1+3]	# pointer to list of ranges
+define	PRL_SZLIST		Memi[$1+4]	# size of list in INTS
+define	PRL_LP			Memi[$1+5]	# offset to next space in list
+
diff --git a/pkg/images/imfit/src/pixlist.x b/pkg/images/imfit/src/pixlist.x
new file mode 100644
index 00000000..066637fd
--- /dev/null
+++ b/pkg/images/imfit/src/pixlist.x
@@ -0,0 +1,369 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "pixlist.h"
+
+.help pixels xtools "Pixel List Handling Tools"
+.nf ________________________________________________________________________
+.fi
+.ih
+PURPOSE
+These routines provide simple pixel list handling facilities and are
+intended as a temporary facility pending full scale completion of
+image masking. The list is stored in the form of ranges as a function
+of line number. Each image line has a offset which may be NULL for
+no entry or an offset into the list itself. The actual list is a set of
+ranges with the ranges for each line delimited by a NULL. Routines
+exist to fetch the ranges for a given line, add or append ranges to a
+given line, fetch the next or previous line number with a non-NULL
+range and specify whether two lines have the same ranges. At present
+the list can grow indefinitely, with additional memory being added as
+necessary. No attempt is made to clean up redundant entries though
+such a faclity could easily be added. The ranges arguments conform
+with the design of the ranges routinesr, with each range consisting
+of and intitial and final entry and a step size. A list of ranges
+is terminated with a NULL
+.ih
+PROCEDURE
+.nf
+prl_init (pl, ncols, nlines)
+
+	pointer	pl		# pointer to list descriptor
+	int	ncols		# number of image columns
+	int	nlines		# number of image lines
+
+nranges = prl_get_ranges (pl, lineno, ranges, max_nranges)
+
+	pointer	pl		# pointer to list descriptor
+	int	lineno		# line number of ranges to be fetched
+	int	ranges[ARB]	# ranges to be output
+	int	max_nranges	# the maximum number of ranges to be output
+
+prl_put_ranges (pl, linemin, linemax, ranges)
+
+	pointer	pl		# pointer to list descriptor
+	int	linemin		# minimum line number
+	int	linemax		# maximum line number
+	int	ranges[ARB]	# ranges to be added to list
+
+prl_append_ranges (pl, linemin, linemax, ranges)
+
+	pointer	pl		# pointer to list descriptor
+	int	linemin		# minimum line number
+	int	linemax		# maximum line number
+	int	ranges[ARB]	# ranges to be added to list
+
+next_lineno/EOF = prl_nextlineno (pl, current_lineno)
+
+	pointer	pl		# pointer to list descriptor
+	int	current_lineno	# current line number
+
+prev_lineno/EOF = prl_prevlineno (pl, current_lineno)
+
+	pointer	pl		# pointer to the list descriptor
+	int	current_lineno	# current line number
+
+YES/NO = prl_eqlines (pl, line1, line2)
+
+	pointer	pl		# pointer to the list descriptor
+	int	line1		# first line number
+	int	line2		# second line number
+
+prl_free	(pl)
+
+	pointer	pl		# pointer to list descriptor
+.fi
+.endhelp ________________________________________________________________
+
+
+# PRL_ADD_RANGES --  Procedure to add the ranges for a given range of
+# line numbers to the pixel list. The new ranges will be appended to any
+# previously existing ranges for the specified line numbers.
+
+procedure prl_add_ranges (pl, linemin, linemax, ranges)
+
+pointer	pl		# pointer to the list descriptor
+int	linemin		# minimum line number
+int	linemax		# maximum line number
+int	ranges[ARB]	# ranges
+
+int	i, j, lc
+int	olp, lp, lnull, lold
+int	nr, nnewr, noldr
+
+begin
+	# check conditions
+	if ((linemin < 1) || (linemax > PRL_NLINES(pl)) || linemin > linemax)
+	    return
+
+	# calculate the length of the range to be appended minus the null
+	nr = 0
+	while (ranges[nr+1] != NULL)
+	    nr = nr + 1
+
+
+	lc = 1
+	olp = -1
+	do i = linemin, linemax {
+
+	    # get offset for line i
+	    lp = Memi[PRL_LINES(pl)+i-1]
+
+	    # if line pointer is undefined
+	    if (lp == NULL) {
+
+		if (lc == 1) {
+
+		    # set line pointer and store
+	    	    Memi[PRL_LINES(pl)+i-1] = PRL_LP(pl)
+		    lnull = PRL_LP(pl)
+
+		    # check the size of the list
+		    if (PRL_SZLIST(pl) < (nr + PRL_LP(pl))) {
+	    		PRL_SZLIST(pl) = PRL_SZLIST(pl) + nr + 1 
+	    		call realloc (PRL_LIST(pl), PRL_SZLIST(pl), TY_INT)
+		    }
+
+		    # move ranges and reset pointers
+		    call amovi (ranges, Memi[PRL_LIST(pl)+PRL_LP(pl)-1], nr)
+		    PRL_LP(pl) = PRL_LP(pl) + nr + 1
+		    Memi[PRL_LIST(pl)+PRL_LP(pl)-2] = NULL
+		    lc = lc + 1
+
+		} else
+
+		    # set line pointer
+		    Memi[PRL_LINES(pl)+i-1] = lnull
+
+	    } else {
+
+		if (lp != olp) {
+
+		    # set line pointer and store
+	    	    Memi[PRL_LINES(pl)+i-1] = PRL_LP(pl)
+		    lold = PRL_LP(pl)
+
+		    # find length of previously defined range and calculate
+		    # length of new ranges
+		    for (j = lp; Memi[PRL_LIST(pl)+j-1] != NULL; j = j + 1)
+		        ;
+		    noldr = j - lp
+		    nnewr = noldr + nr
+
+		    # check size of list
+	            if (PRL_SZLIST(pl) < (nnewr + PRL_LP(pl))) {
+		        PRL_SZLIST(pl) = PRL_SZLIST(pl) + nnewr + 1
+	                call realloc (PRL_LIST(pl), PRL_SZLIST(pl), TY_INT)
+		    }
+
+		    # add ranges to list and update pointers
+		    call amovi (Memi[PRL_LIST(pl)+lp-1],
+		        Memi[PRL_LIST(pl)+PRL_LP(pl)-1], noldr)
+		    PRL_LP(pl) = PRL_LP(pl) + noldr
+		    call amovi (ranges, Memi[PRL_LIST(pl)+PRL_LP(pl)-1], nr)
+		    PRL_LP(pl) = PRL_LP(pl) + nr + 1
+		    Memi[PRL_LIST(pl)+PRL_LP(pl)-2] = NULL
+
+		} else
+
+		    # set line pointers
+		    Memi[PRL_LINES(pl)+i-1] = lold
+
+		olp = lp
+	    }
+	}
+
+end
+
+# PRL_EQLINES -- Routine to test whether two lines have equal ranges.
+# The routine returns YES or NO.
+
+int procedure prl_eqlines (pl, line1, line2)
+
+pointer	pl	# pointer to the list
+int	line1	# line numbers
+int	line2
+
+begin
+	if (Memi[PRL_LINES(pl)+line1-1] == Memi[PRL_LINES(pl)+line2-1])
+	    return (YES)
+	else
+	    return (NO)
+end
+
+# PRL_GET_RANGES -- Procedure to fetch the ranges for the specified lineno.
+# Zero is returned if there are no ranges otherwise the number of ranges
+# are returned. The ranges are stored in an integer array. Three positive
+# numbers are used to define a range a minimum, maximum and a step size.
+# The ranges are delimited by a NULL.
+
+int procedure prl_get_ranges (pl, lineno, ranges, max_nranges)
+
+pointer	pl		# pointer to the pixel list descriptor
+int	lineno		# line number
+int	ranges[ARB]	# array of ranges
+int	max_nranges	# the maximum number of ranges
+
+int	lp, ip
+int	nranges
+
+begin
+	# check for existence of ranges
+	if (Memi[PRL_LINES(pl)+lineno-1] == NULL) {
+	    ranges[1] = NULL
+	    return (0)
+	}
+
+	# set pointer to the first element in list for line lineno
+	lp = PRL_LIST(pl) + Memi[PRL_LINES(pl)+lineno-1] - 1
+
+	# get ranges
+	nranges = 0
+	ip = 1
+	while (Memi[lp+ip-1] != NULL && nranges <= 3 * max_nranges) {
+	    ranges[ip] = Memi[lp+ip-1]
+	    ip = ip + 1
+	    nranges = nranges + 1
+	}
+	ranges[ip] = NULL
+
+	# return nranges
+	if (nranges == 0)
+	    return (nranges)
+	else
+	    return (nranges / 3)
+end
+
+# PRL_NEXTLINENO -- Procedure to fetch the next line number with a set of
+# defined ranges given the current line number. Note that the current
+# line number is updated.
+
+int procedure prl_nextlineno (pl, current_lineno)
+
+pointer	pl			# pointer to the pixel list descriptor
+int	current_lineno		# current line number
+
+int findex, lp
+
+begin
+	findex = max (1, current_lineno + 1)
+	do lp = findex, PRL_NLINES(pl) {
+	    if (Memi[PRL_LINES(pl)+lp-1] != NULL) {
+		current_lineno = lp
+		return (lp)
+	    }
+	}
+
+	return (EOF)
+end
+
+# PRL_PREVLINENO -- Procedure to fetch the first previous line number
+# with a set of defined ranges given the current line number.
+# Note that the current line number is updated.
+
+int procedure prl_prevlineno (pl, current_lineno)
+
+pointer	pl		# pointer to the pixel list descriptor
+int	current_lineno	# current line number
+
+int	findex, lp
+
+begin
+	findex = min (current_lineno - 1, PRL_NLINES(pl))
+	do lp = findex, 1, -1 {
+	    if (Memi[PRL_LINES(pl)+lp-1] != NULL) {
+		current_lineno = lp
+		return (lp)
+	    }
+	}
+
+	return (EOF)
+end
+
+# PRL_PUT_RANGES --  Procedure to add the ranges for a given range of
+# lines to the pixel list. Note that any previously defined ranges are
+# lost.
+
+procedure prl_put_ranges (pl, linemin, linemax, ranges)
+
+pointer	pl		# pointer to the list
+int	linemin		# minimum line
+int	linemax		# maximum line
+int	ranges[ARB]	# list of ranges
+
+int	i
+int	len_range
+
+begin
+	# check boundary conditions
+	if ((linemin < 1) || (linemax > PRL_NLINES(pl)) || (linemin > linemax))
+	    return
+
+	# determine length of range string minus the NULL
+	len_range = 0
+	while (ranges[len_range+1] != NULL)
+	    len_range = len_range + 1
+
+	# check space allocation
+	if (PRL_SZLIST(pl) < (len_range + PRL_LP(pl))) {
+	    PRL_SZLIST(pl) = PRL_SZLIST(pl) + len_range + 1
+	    call realloc (PRL_LIST(pl), PRL_SZLIST(pl), TY_INT)
+	}
+
+	# set the line pointers
+	do i = linemin, linemax
+	    Memi[PRL_LINES(pl)+i-1] = PRL_LP(pl)
+
+	# add ranges
+	call amovi (ranges, Memi[PRL_LIST(pl)+PRL_LP(pl)-1], len_range)
+	PRL_LP(pl) = PRL_LP(pl) + len_range + 1
+	Memi[PRL_LIST(pl)+PRL_LP(pl)-2] = NULL
+end
+
+
+# PLR_FREE -- Procedure to free the pixel list descriptor
+
+procedure prl_free (pl)
+
+pointer	pl		# pointer to pixel list descriptor
+
+begin
+	if (pl == NULL)
+	    return
+
+	if (PRL_LIST(pl) != NULL)
+	    call mfree (PRL_LIST(pl), TY_INT)
+	if (PRL_LINES(pl) != NULL)
+	    call mfree (PRL_LINES(pl), TY_INT)
+
+	call mfree (pl, TY_STRUCT)
+end
+
+# PRL_INIT -- Procedure to initialize the pixel list. Ncols and nlines are
+# the number of columns and lines respectively in the associated IRAF
+# image.
+
+procedure prl_init (pl, ncols, nlines)
+
+pointer	pl		# pixel list descriptor
+int	ncols		# number of image columns 
+int	nlines		# number of image lines
+
+begin
+	# allocate space for a pixel list descriptor
+	call malloc (pl, LEN_PLSTRUCT, TY_STRUCT)
+
+	# initialize
+	PRL_NCOLS(pl) = ncols
+	PRL_NLINES(pl) = nlines
+
+	# allocate space for the line pointers
+	call malloc (PRL_LINES(pl), PRL_NLINES(pl), TY_INT) 
+	call amovki (NULL, Memi[PRL_LINES(pl)], PRL_NLINES(pl))
+
+	# set pointer to next free element
+	PRL_LP(pl) = 1
+
+	# allocate space for the actual list
+	call malloc (PRL_LIST(pl), PRL_NLINES(pl), TY_INT)
+	PRL_SZLIST(pl) = PRL_NLINES(pl)
+end
diff --git a/pkg/images/imfit/src/ranges.x b/pkg/images/imfit/src/ranges.x
new file mode 100644
index 00000000..19dc5c0e
--- /dev/null
+++ b/pkg/images/imfit/src/ranges.x
@@ -0,0 +1,524 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	<ctype.h>
+
+.help ranges xtools "Range Parsing Tools"
+.ih
+PURPOSE
+
+These tools
+parse a string using a syntax to represent integer values, ranges, and
+steps.   The parsed string is used to generate a list of integers for various
+purposes such as specifying lines or columns in an image or tape file numbers.
+.ih
+SYNTAX
+
+The syntax for the range string consists of positive integers, '-' (minus),
+'x', ',' (comma), and whitespace.  The commas and whitespace are ignored
+and may be freely used for clarity.  The remainder of the string consists
+of sequences of five fields.  The first field is the beginning of a range,
+the second is a '-', the third is the end of the range, the fourth is
+a 'x', and the fifth is a step size.  Any of the five fields may be
+missing causing various default actions.  The defaults are illustrated in
+the following table.
+
+.nf
+-3x1	A missing starting value defaults to 1.
+2-x1	A missing ending value defaults to MAX_INT.
+2x1	A missing ending value defaults to MAX_INT.
+2-4	A missing step defaults to 1.
+4	A missing ending value and step defaults to an ending
+	value equal to the starting value and a step of 1.
+x2	Missing starting and ending values defaults to
+	the range 1 to MAX_INT with the specified step.
+""	The null string is equivalent to "1 - MAX_INT x 1",
+	i.e all positive integers.
+.fi
+
+The specification of several ranges yields the union of the ranges.
+.ih
+EXAMPLES
+
+The following examples further illustrate the range syntax.
+
+.nf
+-	All positive integers.
+1,5,9	A list of integers equivalent to 1-1x1,5-5x1,9-9x1.
+x2	Every second positive integer starting with 1.
+2x3	Every third positive integer starting with 2.
+-10	All integers between 1 and 10.
+5-	All integers greater than or equal to 5.
+9-3x1	The integers 3,6,9.
+.fi
+.ih
+PROCEDURES
+
+.ls 4 is_decode_ranges
+
+.nf
+int procedure is_decode_ranges (range_string, ranges, max_ranges, minimum,
+    maximum, nvalues)
+
+char	range_string[ARB]	# Range string to be decoded
+int	ranges[3, max_ranges]	# Range array
+int	max_ranges		# Maximum number of ranges
+int	minimum, maximum	# Minimum and maximum range values allowed
+int	nvalues			# The number of values in the ranges
+.fi
+
+The range string is decoded into an integer array of maximum dimension
+3 * max_ranges.  Each range consists of three consecutive integers
+corresponding to the starting and ending points of the range and the
+step size.  The number of integers covered by the ranges is returned
+as nvalue.  The end of the set of ranges is marked by a NULL.
+The returned status is either ERR or OK.
+.le
+.ls 4 is_next_number, is_previous_number
+
+.nf
+int procedure is_next_number (ranges, number)
+int procedure is_previous_number (ranges, number)
+
+int	ranges[ARB]		# Range array
+int	number			# Both input and output parameter
+.fi
+
+Given a value for number the procedures find the next (previous) number in
+increasing (decreasing)
+value within the set of ranges.  The next (previous) number is returned in
+the number argument.  A returned status is either OK or EOF.
+EOF indicates that there are no greater values.  The usual usage would
+be in a loop of the form:
+
+.nf
+	number = 0
+	while (is_next_number (ranges, number) != EOF) {
+	    <Statements using number>
+	}
+.fi
+.le
+.ls 4 is_in_rangelist
+
+.nf
+bool procedure is_in_rangelist (ranges, number)
+
+int	ranges[ARB]		# Ranges array
+int	number			# Number to check againts ranges
+.fi
+
+A boolean value is returned indicating whether number is covered by
+the ranges.
+
+.endhelp
+
+
+# IS_DECODE_RANGES -- Parse a string containing a list of integer numbers or
+# ranges, delimited by either spaces or commas.  Return as output a list
+# of ranges defining a list of numbers, and the count of list numbers.
+# Range limits must be positive nonnegative integers.  ERR is returned as
+# the function value if a conversion error occurs.  The list of ranges is
+# delimited by a single NULL.
+
+
+int procedure is_decode_ranges (range_string, ranges, max_ranges, minimum,
+    maximum, nvalues)
+
+char	range_string[ARB]	# Range string to be decoded
+int	ranges[3, max_ranges]	# Range array
+int	max_ranges		# Maximum number of ranges
+int	minimum, maximum	# Minimum and maximum range values allowed
+int	nvalues			# The number of values in the ranges
+
+int	ip, nrange, out_of_range, a, b, first, last, step, ctoi()
+
+begin
+	ip = 1
+	nrange = 1
+	nvalues = 0
+	out_of_range = 0
+
+	while (nrange < max_ranges) {
+	    # Default values
+	    a = minimum
+	    b = maximum
+	    step = 1
+
+	    # Skip delimiters
+	    while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+		ip = ip + 1
+
+	    # Get first limit.
+	    # Must be a number, '*', '-', 'x', or EOS.  If not return ERR.
+	    if (range_string[ip] == EOS) {			# end of list
+		if (nrange == 1) {
+		    if (out_of_range == 0) {
+		        # Null string defaults
+		        ranges[1, 1] = a
+		        ranges[2, 1] = b
+		        ranges[3, 1] = step
+		        ranges[1, 2] = NULL
+	    	        nvalues = (b - a) / step + 1
+		        return (OK)
+		    } else {
+			# Only out of range data
+			return (ERR)
+		    }
+		} else {
+		    ranges[1, nrange] = NULL
+		    return (OK)
+		}
+	    } else if (range_string[ip] == '-')
+		;
+	    else if (range_string[ip] == '*')
+		;
+	    else if (range_string[ip] == 'x')
+		;
+	    else if (IS_DIGIT(range_string[ip])) {		# ,n..
+		if (ctoi (range_string, ip, a) == 0)
+		    return (ERR)
+	    } else
+		return (ERR)
+
+	    # Skip delimiters
+	    while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+		ip = ip + 1
+
+	    # Get last limit
+	    # Must be '-', '*', or 'x' otherwise b = a.
+	    if (range_string[ip] == 'x')
+		;
+	    else if ((range_string[ip] == '-') || (range_string[ip] == '*')) {
+		ip = ip + 1
+	        while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+		    ip = ip + 1
+		if (range_string[ip] == EOS)
+		    ;
+		else if (IS_DIGIT(range_string[ip])) {
+		    if (ctoi (range_string, ip, b) == 0)
+		        return (ERR)
+		} else if (range_string[ip] == 'x')
+		    ;
+		else
+		    return (ERR)
+	    } else
+		b = a
+
+	    # Skip delimiters
+	    while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+		ip = ip + 1
+
+	    # Get step.
+	    # Must be 'x' or assume default step.
+	    if (range_string[ip] == 'x') {
+		ip = ip + 1
+	        while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+		    ip = ip + 1
+		if (range_string[ip] == EOS)
+		    ;
+		else if (IS_DIGIT(range_string[ip])) {
+		    if (ctoi (range_string, ip, step) == 0)
+		        ;
+		} else if (range_string[ip] == '-')
+		    ;
+		else if (range_string[ip] == '*')
+		    ;
+		else
+		    return (ERR)
+	    }
+
+	    # Output the range triple.
+	    first = min (a, b)
+	    last = max (a, b)
+	    if (first < minimum)
+		first = minimum + mod (step - mod (minimum - first, step), step)
+	    if (last > maximum)
+		last = maximum - mod (last - maximum, step)
+	    if (first <= last) {
+	        ranges[1, nrange] = first
+	        ranges[2, nrange] = last
+	        ranges[3, nrange] = step
+	        nvalues = nvalues + (last - first) / step + 1
+		nrange = nrange + 1
+	    } else
+		out_of_range = out_of_range + 1
+	}
+
+	return (ERR)					# ran out of space
+end
+
+
+# IS_NEXT_NUMBER -- Given a list of ranges and the current file number,
+# find and return the next file number.  Selection is done in such a way
+# that list numbers are always returned in monotonically increasing order,
+# regardless of the order in which the ranges are given.  Duplicate entries
+# are ignored.  EOF is returned at the end of the list.
+
+int procedure is_next_number (ranges, number)
+
+int	ranges[ARB]		# Range array
+int	number			# Both input and output parameter
+
+int	ip, first, last, step, next_number, remainder
+
+begin
+	# If number+1 is anywhere in the list, that is the next number,
+	# otherwise the next number is the smallest number in the list which
+	# is greater than number+1.
+
+	number = number + 1
+	next_number = MAX_INT
+
+	for (ip=1;  ranges[ip] != NULL;  ip=ip+3) {
+	    first = ranges[ip]
+	    last = ranges[ip+1]
+	    step = ranges[ip+2]
+	    if (number >= first && number <= last) {
+		remainder = mod (number - first, step)
+		if (remainder == 0)
+		    return (number)
+		if (number - remainder + step <= last)
+		    next_number = number - remainder + step
+	    } else if (first > number)
+		next_number = min (next_number, first)
+	}
+
+	if (next_number == MAX_INT)
+	    return (EOF)
+	else {
+	    number = next_number
+	    return (number)
+	}
+end
+
+
+# IS_PREVIOUS_NUMBER -- Given a list of ranges and the current file number,
+# find and return the previous file number.  Selection is done in such a way
+# that list numbers are always returned in monotonically decreasing order,
+# regardless of the order in which the ranges are given.  Duplicate entries
+# are ignored.  EOF is returned at the end of the list.
+
+int procedure is_previous_number (ranges, number)
+
+int	ranges[ARB]		# Range array
+int	number			# Both input and output parameter
+
+int	ip, first, last, step, next_number, remainder
+
+begin
+	# If number-1 is anywhere in the list, that is the previous number,
+	# otherwise the previous number is the largest number in the list which
+	# is less than number-1.
+
+	number = number - 1
+	next_number = 0
+
+	for (ip=1;  ranges[ip] != NULL;  ip=ip+3) {
+	    first = ranges[ip]
+	    last = ranges[ip+1]
+	    step = ranges[ip+2]
+	    if (number >= first && number <= last) {
+		remainder = mod (number - first, step)
+		if (remainder == 0)
+		    return (number)
+		if (number - remainder >= first)
+		    next_number = number - remainder
+	    } else if (last < number) {
+		remainder = mod (last - first, step)
+		if (remainder == 0)
+		    next_number = max (next_number, last)
+		else if (last - remainder >= first)
+		    next_number = max (next_number, last - remainder)
+	    }
+	}
+
+	if (next_number == 0)
+	    return (EOF)
+	else {
+	    number = next_number
+	    return (number)
+	}
+end
+
+
+# IS_IN_RANGELLIST -- Test number to see if it is in range.
+
+bool procedure is_in_rangelist (ranges, number)
+
+int	ranges[ARB]		# Range array
+int	number			# Number to be tested against ranges
+
+int	ip, first, last, step
+
+begin
+	for (ip=1;  ranges[ip] != NULL;  ip=ip+3) {
+	    first = ranges[ip]
+	    last = ranges[ip+1]
+	    step = ranges[ip+2]
+	    if (number >= first && number <= last)
+		if (mod (number - first, step) == 0)
+		    return (TRUE)
+	}
+
+	return (FALSE)
+end
+
+
+# IS_EXPAND_RANGES -- Expand a range string into a array of values.
+
+int procedure is_expand_ranges (ranges, array, max_nvalues)
+
+int	ranges[ARB]			# Range array
+int	array[max_nvalues]		# Array of values
+int	max_nvalues			# Maximum number of values
+
+int	n, value
+
+int	is_next_number()
+
+begin
+	n = 0
+	value = 0
+	while ((n < max_nvalues) && (is_next_number (ranges, value) != EOF)) {
+	    n = n + 1
+	    array[n] = value
+	}
+
+	return (n)
+end
+
+
+# IS_SELECT_RANGES -- Select array values in the ranges.
+# The input and output arrays may be the same.
+
+procedure is_select_ranges (a, b, ranges)
+
+real	a[ARB]				# Input array
+real	b[ARB]				# Output array
+int	ranges[3, ARB]			# Ranges
+
+int	i, j, npts, nmove
+
+begin
+	npts = 0
+	for (i = 1; ranges[1, i] != NULL; i = i + 1) {
+	    if (ranges[3, i] == 1) {
+	        nmove = ranges[2, i] - ranges[1, i] + 1
+	        call amovr (a[ranges[1, i]], b[npts + 1], nmove)
+		npts = npts + nmove
+	    } else {
+		do j = ranges[1, i], ranges[2, i], ranges[3, i] {
+		    npts = npts + 1
+		    b[npts] = a[j]
+		}
+	    }
+	}
+end
+
+
+# IS_CHOOSE_RANGESI -- Copy the selected values from array a to b.
+
+int procedure is_choose_rangesi (indices, a, b, npts, ifirst, ilast)
+
+int	indices[ARB]		# array of indices
+int	a[ARB]			# input array
+int	b[ARB]			# output array
+int	npts			# number of points
+int	ifirst			# first index
+int	ilast			# last index
+
+int	i, element
+
+begin
+	element = 1
+	do i = 1, npts {
+	    if (indices[i] < ifirst || indices[i] > ilast)
+		next
+	    b[element] = a[indices[i]]
+	    element = element + 1
+	}
+	return (element - 1)
+end
+
+
+# IS_CHOOSE_RANGESR -- Copy the selected values from array a to b.
+
+int procedure is_choose_rangesr (indices, a, b, npts, ifirst, ilast)
+
+int	indices[ARB]		# array of indices
+real	a[ARB]			# input array
+real	b[ARB]			# output array
+int	npts			# number of points
+int	ifirst			# first element to be extracted
+int	ilast			# last element to be extracted
+
+int	i, element
+
+begin
+	element = 1
+	do i = 1, npts {
+	    if (indices[i] < ifirst || indices[i] > ilast)
+		next
+	    b[element] = a[indices[i]]
+	    element = element + 1
+	}
+	return (element - 1)
+end
+
+
+# IS_MAKE_RANGES -- Procedure to make a set of ranges from an ordered list
+# of column numbers. Only a step size of 1 is checked for.
+
+int procedure is_make_ranges (list, npts, ranges, max_nranges)
+
+int	list[ARB]	# list of column numbers in increasing order
+int	npts		# number of list elements
+int	ranges[ARB]	# output ranges
+int	max_nranges	# the maximum number of ranges
+
+bool	next_range
+int	ip, op, nranges
+
+begin
+	# If zero list elements return
+	if (npts == 0) {
+	    ranges[1] = NULL	
+	    return (0)
+	}
+
+	# Initialize
+	nranges = 0
+	ranges[1] = list[1]
+	op = 2
+	next_range = false
+
+	# Loop over column list
+	for (ip = 2; ip <= npts && nranges < max_nranges; ip = ip + 1) {
+	    if ((list[ip] != (list[ip-1] + 1))) {
+		ranges[op] = list[ip-1]
+		op = op + 1
+		ranges[op] = 1
+		op = op + 1
+		nranges = nranges + 1
+		ranges[op] = list[ip]
+		op = op + 1
+	    }
+	}
+
+	# finish off
+	if (npts == 1) {
+	    ranges[op] = list[npts]
+	    ranges[op+1] = 1
+	    ranges[op+2] = NULL
+	    nranges = 1
+	} else if (nranges == max_nranges) {
+	    ranges[op-1] = NULL
+	} else {
+	    ranges[op] = list[npts]
+	    ranges[op+1] = 1
+	    ranges[op+2] = NULL
+	    nranges = nranges + 1
+	}
+
+	return (nranges)
+end
diff --git a/pkg/images/imfit/src/t_imsurfit.x b/pkg/images/imfit/src/t_imsurfit.x
new file mode 100644
index 00000000..2a93b8b2
--- /dev/null
+++ b/pkg/images/imfit/src/t_imsurfit.x
@@ -0,0 +1,400 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <error.h>
+include <imhdr.h>
+include "imsurfit.h"
+
+# T_IMSURFIT -- Fit a surface function to an image
+#
+# 1. A user selected function is fit to each surface.
+# 2. Only the selected regions of the image are fit.
+# 3. Deviant pixels may be rejected from the fit.
+# 4. The user selects the type of output image. The choices are:
+#    a. the fitted image.
+#    b. the input image with deviant  pixels replaced by
+#       the fitted values
+#    c. the input image minus the fitted image.
+#    d. the ratio of the input image and the fit where
+#       pixels less than div_min are set to a ratio of 1.
+
+
+procedure t_imsurfit ()
+
+char	imtlist1[SZ_LINE]			# Input image list
+char	imtlist2[SZ_LINE]			# Output image list
+char	image1[SZ_FNAME]			# Input image
+char	image2[SZ_FNAME]			# Output image
+
+char	str[SZ_LINE], region_str[SZ_LINE]
+int	list1, list2, region_type
+pointer	im1, im2, imfit, gl, sp
+
+bool	clgetb()
+int	imtopen(), imtgetim(), imtlen(), btoi(), clgeti(), clgwrd()
+pointer	immap()
+real	clgetr()
+
+begin
+	# Allocate space for imfit structure.
+	call smark (sp)
+	call salloc (imfit, LEN_IMSFSTRUCT, TY_STRUCT)
+
+	# Get task parameters.
+	call clgstr ("input", imtlist1, SZ_FNAME)
+	call clgstr ("output", imtlist2, SZ_FNAME)
+	TYPE_OUTPUT(imfit) = clgwrd ("type_output", str, SZ_LINE,
+		      	     ",fit,clean,residual,response,")
+	DIV_MIN(imfit) = clgetr ("div_min")
+
+	# Get surface ftting parameters.
+	SURFACE_TYPE(imfit) = clgwrd ("function", str, SZ_LINE,
+			      ",legendre,chebyshev,spline3,spline1,")
+	XORDER(imfit) = clgeti ("xorder")
+	YORDER(imfit) = clgeti ("yorder")
+	CROSS_TERMS(imfit) = btoi (clgetb ("cross_terms"))
+
+	# Get median processing parameters.
+	XMEDIAN(imfit) = clgeti ("xmedian")
+	YMEDIAN(imfit) = clgeti ("ymedian")
+	MEDIAN_PERCENT(imfit) = clgetr ("median_percent")
+	if (XMEDIAN(imfit) > 1 || YMEDIAN(imfit) > 1)
+	    MEDIAN(imfit) = YES
+	else
+	    MEDIAN(imfit) = NO
+
+	# Get rejection cycle parameters.
+	NITER(imfit) = clgeti ("niter")
+	LOWER(imfit) = clgetr ("lower")
+	UPPER(imfit) = clgetr ("upper")
+	NGROW(imfit) = clgeti ("ngrow")
+
+	if (MEDIAN(IMFIT) == YES) {
+	    REJECT(imfit) = NO
+	    NITER(imfit) = 0
+	} else if (NITER(imfit) > 0 && (LOWER(imfit) > 0. || UPPER(imfit) > 0.))
+	    REJECT(imfit) = YES
+	else {
+	    REJECT(imfit) = NO
+	    NITER(imfit) = 0
+	}
+
+	# Checking sigmas for cleaning.
+	if (TYPE_OUTPUT(imfit) == CLEAN && MEDIAN(imfit) == YES)
+	    call error (0,
+	    "T_IMSURFIT: Clean option and median processing are exclusive.")
+	if (TYPE_OUTPUT(imfit) == CLEAN && NITER(imfit) <= 0)
+	    call error (0, "T_IMSURFIT: Clean option requires non-zero niter.")
+	if (TYPE_OUTPUT(imfit) == CLEAN && LOWER(imfit) <= 0. &&
+	    UPPER(imfit) <= 0.)
+	    call error (0, "T_IMSURFIT: Clean option requires non-zero sigma.")
+
+	# Get regions to be fit.
+	gl = NULL
+	region_type = clgwrd ("regions", str, SZ_LINE,
+		      ",all,columns,rows,border,sections,circle,invcircle,")
+	switch (region_type) {
+	case ALL:
+	    ;
+	case BORDER:
+	    call clgstr ("border", region_str, SZ_LINE)
+	case SECTIONS:
+	    call clgstr ("sections", region_str, SZ_LINE)
+	case COLUMNS:
+	    call clgstr ("columns", region_str, SZ_LINE)
+	case ROWS:
+	    call clgstr ("rows", region_str, SZ_LINE)
+	case CIRCLE, INVCIRCLE:
+	    call clgstr ("circle", region_str, SZ_LINE)
+	}
+
+	# Expand the input and output image lists.
+	list1 = imtopen (imtlist1)
+	list2 = imtopen (imtlist2)
+	if (imtlen (list1) != imtlen (list2)) {
+	    call imtclose (list1)
+	    call imtclose (list2)
+	    call error (0, "Number of input and output images not the same.")
+	}
+
+	# Do each set of input and output images.
+	while ((imtgetim (list1, image1, SZ_FNAME) != EOF) &&
+	      (imtgetim (list2, image2, SZ_FNAME) != EOF)) {
+	    
+	    im1 = immap (image1, READ_ONLY, 0)
+	    im2 = immap (image2, NEW_COPY, im1)
+
+	    iferr {
+	        if (region_type != ALL)
+		    call make_good_list (im1, gl, region_type, region_str)
+	        call imsurfit (im1, im2, imfit, gl)
+	    } then
+		call erract (EA_WARN)
+
+	    call imunmap (im1)
+	    call imunmap (im2)
+	    call prl_free (gl)
+	}
+
+	# Cleanup.
+	call sfree (sp)
+	call imtclose (list1)
+	call imtclose (list2)
+end
+
+
+# MAKE_GOOD_LIST -- Procedure to make a list of good regions. The program
+# returns an error message if no good regions are defined. The good
+# list parameter is set to NULL if the whole image is to be fit. This routine
+# uses both the ranges and pixlist package which will be replaced by image
+# masking.
+
+procedure make_good_list (im, gl, region_type, region_string)
+
+pointer	im			# pointer to the image
+pointer	gl			# good pixel list descriptor
+int	region_type		# type of good region
+char	region_string[ARB]	# region parameters
+
+int	i, ip, zero, nvals, range_min, r2, xdist, max_nranges
+int	x1, x2, y1, y2, temp, border, xcenter, ycenter, radius
+int	columns[7]
+pointer	sp, ranges, list
+
+bool	is_in_rangelist()
+int	is_next_number(), is_decode_ranges(), open(), fscan(), nscan(), ctoi()
+errchk	open, close
+
+begin
+	# Determine the maximum number of images.
+	max_nranges = IM_LEN(im,1)
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (ranges, 3 * max_nranges + 1, TY_INT)
+
+	# Compute the good pixel list.
+	switch (region_type) {
+	case ROWS:
+
+	    # Decode the row ranges.
+	    if (is_decode_ranges (region_string, Memi[ranges], max_nranges, 1,
+	        int (IM_LEN(im,2)), nvals) == ERR)
+		call error (0, "MAKE_GOOD_LIST: Error decoding row string.")
+	    if (nvals == 0) 
+		call error (0, "MAKE_GOOD_LIST: no good rows.") 
+	    if (nvals == IM_LEN(im,2)) {
+		call sfree (sp)
+		return
+	    }
+
+	    # Intialize the good pixel list.
+	    call prl_init (gl, int (IM_LEN(im,1)), int (IM_LEN(im,2)))
+
+	    # Set column limits using the ranges format.
+	    columns[1] = 1
+	    columns[2] = IM_LEN(im,1)
+	    columns[3] = 1
+	    columns[4] = NULL
+
+	    # Set column limits for the specied lines.
+	    zero = 0
+	    range_min = is_next_number (Memi[ranges], zero)
+	    while (range_min != EOF) {
+	        for (i = range_min; i <= IM_LEN(im,2) + 1; i = i + 1) {
+		    if (!is_in_rangelist (Memi[ranges], i) ||
+		        i == IM_LEN(im,2)+1) {
+		        call prl_put_ranges (gl, range_min, i-1, columns)
+		        break
+		    }
+		}
+		range_min = is_next_number (Memi[ranges], i)
+	    }
+
+	case COLUMNS:
+
+	    # Set the specified columns.
+	    if (is_decode_ranges (region_string, Memi[ranges], max_nranges, 1,
+	        int (IM_LEN(im,1)), nvals) == ERR)
+		call error (0, "MAKE_GOOD_LIST: Error decoding column string.")
+	    if (nvals == 0)
+		call error (0, "MAKE_GOOD_LIST: No good columns.")
+	    if (nvals == IM_LEN(im,1)) {
+		call sfree (sp)
+		return
+	    }
+
+	    # Make the good pixel list.
+	    call prl_init (gl, int (IM_LEN(im,1)), int (IM_LEN(im,2)))
+	    call prl_add_ranges (gl, 1, int (IM_LEN(im,2)), Memi[ranges])
+
+	case CIRCLE, INVCIRCLE:
+	    
+	    # Get the parameters of the circle.
+	    ip = 1
+	    if (ctoi (region_string, ip, xcenter) <= 0)
+		call error (0, "MAKE_GOOD_LIST: Error decoding xcenter.")
+	    if (ctoi (region_string, ip, ycenter) <= 0)
+		call error (0, "MAKE_GOOD_LIST: Error decoding ycenter.")
+	    if (ctoi (region_string, ip, radius) <= 0)
+		call error (0, "MAKE_GOOD_LIST: Error decoding radius.")
+
+	    y1 = max (1, ycenter - radius)
+	    y2 = min (int (IM_LEN(im,2)), ycenter + radius)
+	    x1 = max (1, xcenter - radius)
+	    x2 = min (int (IM_LEN(im,1)), xcenter + radius)
+	    if (region_type == CIRCLE) {
+	        if (y1 > IM_LEN(im,2) || y2 < 1 || x1 > IM_LEN(im,1) || x2 < 1)
+		    call error (0, "MAKE_GOOD_LIST: No good regions.")
+	    }
+
+	    # Create the good pixel list.
+	    call prl_init (gl, int (IM_LEN(im,1)), int (IM_LEN(im,2)))
+
+	    r2 = radius ** 2
+	    if (region_type == CIRCLE) {
+	        do i = y1, y2 {
+		    xdist = sqrt (real (r2 - (ycenter - i) ** 2))
+		    x1 = max (1, xcenter - xdist)
+		    x2 = min (IM_LEN(im,1), xcenter + xdist)
+		    columns[1] = x1
+		    columns[2] = x2
+		    columns[3] = 1
+		    columns[4] = NULL
+		    call prl_put_ranges (gl, i, i, columns)
+	        }
+	    } else if (region_type == INVCIRCLE) {
+		do i = 1, y1 - 1 {
+		    columns[1] = 1
+		    columns[2] = IM_LEN(im,1)
+		    columns[3] = 1
+		    columns[4] = NULL
+		    call prl_put_ranges (gl, i, i, columns)
+		}
+		do i = y2 + 1, IM_LEN(im,2) {
+		    columns[1] = 1
+		    columns[2] = IM_LEN(im,1)
+		    columns[3] = 1
+		    columns[4] = NULL
+		    call prl_put_ranges (gl, i, i, columns)
+		}
+		do i = y1, y2 {
+		    xdist = sqrt (real (r2 - (ycenter - i) ** 2))
+		    x1 = max (1, xcenter - xdist)
+		    x2 = min (IM_LEN(im,1), xcenter + xdist)
+		    if (x1 > 1) {
+			columns[1] = 1
+			columns[2] = x1 - 1
+			columns[3] = 1
+			if (x2 < IM_LEN(im,1)) {
+			    columns[4] = x2 + 1
+			    columns[5] = IM_LEN(im,1)
+			    columns[6] = 1
+			    columns[7] = NULL
+			} else
+			    columns[4] = NULL
+		    } else if (x2 < IM_LEN(im,1)) {
+			columns[1] = x2 + 1
+			columns[2] = IM_LEN(im,1)
+			columns[3] = 1
+			columns[4] = NULL
+		    } else
+			columns[1] = NULL
+		    call prl_put_ranges (gl, i, i, columns)
+		}
+	    }
+
+
+	case SECTIONS:
+
+	    # Open file of sections.
+	    list = open (region_string, READ_ONLY, TEXT_FILE)
+	    call prl_init (gl, int (IM_LEN(im,1)), int (IM_LEN(im,2)))
+
+	    # Scan the list.
+	    while (fscan (list) != EOF) {
+		
+		# Fetch parameters from list.
+		call gargi (x1)
+		call gargi (x2)
+		call gargi (y1)
+		call gargi (y2)
+		if (nscan() != 4)
+		    next
+		
+		# Check and correct for out of bounds limits.
+		x1 = max (1, min (IM_LEN(im,1), x1))
+		x2 = min (IM_LEN(im,1), max (1, x2))
+		y1 = max (1, min (IM_LEN(im,2), y1))
+		y2 = min (IM_LEN(im,2), max (1, y2))
+
+		# Check the order.
+		if (x2 < x1) {
+		    temp = x1
+		    x1 = x2
+		    x2 = temp
+		}
+		if (y2 < y1) {
+		    temp = y1
+		    y1 = y2
+		    y2 = temp
+		}
+
+		# If entire image return.
+		if ((x1 == 1) && (x2 == IM_LEN(im,1)) && (y1 == 1) &&
+		    (y2 == IM_LEN(im,2))) {
+		    call prl_free (gl)
+		    gl = NULL
+		    break
+		}
+
+		# Set ranges.
+		columns[1] = x1
+		columns[2] = x2
+		columns[3] = 1
+		columns[4] = NULL
+		call prl_add_ranges (gl, y1, y2, columns)
+	    }
+
+	    call close (list)
+
+	case BORDER:
+
+	    # Decode border parameter.
+	    ip = 1
+	    if (ctoi (region_string, ip, border) == ERR)
+		call error (0, "MAKE_GOOD_LIST: Error decoding border string.")
+	    if (border < 1)
+		call error (0, "MAKE_GOOD_LIST: No border.")
+	    if ((border > IM_LEN(im,1)/2) && (border > IM_LEN(im,2)/2)) {
+		call sfree (sp)
+		return
+	    }
+
+	    # Intialize list.
+	    call prl_init (gl, int (IM_LEN(im,1)), int (IM_LEN(im,2)))
+	    y1 = 1 + border - 1
+	    y2 = IM_LEN(im,2) - border + 1
+	    columns[1] = 1
+	    columns[2] = IM_LEN(im,1)
+	    columns[3] = 1
+	    columns[4] = NULL
+
+	    # Set ranges for top and bottom edges of image.
+	    call prl_put_ranges (gl, 1, y1, columns)
+	    call prl_put_ranges (gl, y2, int (IM_LEN(im,2)), columns)
+
+	    columns[1] = 1
+	    columns[2] = y1
+	    columns[3] = 1
+	    columns[4] = NULL
+	    call prl_put_ranges (gl, y1 + 1, y2 - 1, columns)
+
+	    columns[1] = IM_LEN(im,1) - border + 1
+	    columns[2] = IM_LEN(im,1)
+	    columns[3] = 1
+	    columns[4] = NULL
+	    call prl_add_ranges (gl, y1 + 1, y2 - 1, columns)
+	}
+
+	call sfree (sp)
+end
diff --git a/pkg/images/imfit/src/t_lineclean.x b/pkg/images/imfit/src/t_lineclean.x
new file mode 100644
index 00000000..4acb9752
--- /dev/null
+++ b/pkg/images/imfit/src/t_lineclean.x
@@ -0,0 +1,270 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<pkg/gtools.h>
+
+# LINECLEAN -- Fit a function to the image lines and output an image
+# with rejected points replaced by the fit.  The fitting parameters may be
+# set interactively using the icfit package.
+
+procedure t_lineclean ()
+
+int	listin				# Input image list
+int	listout				# Output image list
+char	sample[SZ_LINE]			# Sample ranges
+int	naverage			# Sample averaging size
+char	function[SZ_LINE]		# Curve fitting function
+int	order				# Order of curve fitting function
+real	low_reject, high_reject		# Rejection threshold
+int	niterate			# Number of rejection iterations
+real	grow				# Rejection growing radius
+bool	interactive			# Interactive?
+
+char	input[SZ_LINE]			# Input image
+char	output[SZ_FNAME]		# Output image
+pointer	in, out				# IMIO pointers
+pointer	ic				# ICFIT pointer
+pointer	gt				# GTOOLS pointer
+
+int	imtopen(), imtgetim(), imtlen(), gt_init()
+int	clgeti()
+real	clgetr()
+bool	clgetb()
+
+begin
+	# Get input and output lists and check that the number of images
+	# are the same.
+
+	call clgstr ("input", input, SZ_LINE)
+	listin = imtopen (input)
+	call clgstr ("output", input, SZ_LINE)
+	listout = imtopen (input)
+	if (imtlen (listin) != imtlen (listout)) {
+	    call imtclose (listin)
+	    call imtclose (listout)
+	    call error (0, "Input and output image lists do not match")
+	}
+
+	# Get task parameters.
+
+	call clgstr ("sample", sample, SZ_LINE)
+	naverage = clgeti ("naverage")
+	call clgstr ("function", function, SZ_LINE)
+	order = clgeti ("order")
+	low_reject = clgetr ("low_reject")
+	high_reject = clgetr ("high_reject")
+	niterate = clgeti ("niterate")
+	grow = clgetr ("grow")
+	interactive = clgetb ("interactive")
+
+
+	# Set the ICFIT pointer structure.
+	call ic_open (ic)
+	call ic_pstr (ic, "sample", sample)
+	call ic_puti (ic, "naverage", naverage)
+	call ic_pstr (ic, "function", function)
+	call ic_puti (ic, "order", order)
+	call ic_putr (ic, "low", low_reject)
+	call ic_putr (ic, "high", high_reject)
+	call ic_puti (ic, "niterate", niterate)
+	call ic_putr (ic, "grow", grow)
+	call ic_pstr (ic, "ylabel", "")
+
+	gt = gt_init()
+	call gt_sets (gt, GTTYPE, "line")
+
+	# Clean the lines in each input image.
+
+	while ((imtgetim (listin, input, SZ_FNAME) != EOF) &&
+	    (imtgetim (listout, output, SZ_FNAME) != EOF)) {
+
+	    call lc_immap (input, output, in, out)
+	    call lineclean (in, out, ic, gt, input, interactive)
+	    call imunmap (in)
+	    call imunmap (out)
+	}
+
+	call ic_closer (ic)
+	call gt_free (gt)
+	call imtclose (listin)
+	call imtclose (listout)
+end
+
+
+# LINECLEAN -- Given the image descriptor determine the fitting function
+# for each line and create an output image.  If the interactive flag
+# is set then set the fitting parameters interactively.
+
+procedure lineclean (in, out, ic, gt, title, interactive)
+
+pointer	in				# IMIO pointer for input image
+pointer	out				# IMIO pointer for output image
+pointer	ic				# ICFIT pointer
+pointer	gt				# GTIO pointer
+char	title[ARB]			# Title
+bool	interactive			# Interactive?
+
+char	graphics[SZ_FNAME]
+int	i, nx, new
+long	inline[IM_MAXDIM], outline[IM_MAXDIM]
+pointer	cv, gp, sp, x, wts, indata, outdata
+
+int	lf_getline(), imgnlr(), impnlr(), strlen()
+pointer	gopen()
+
+begin
+	# Allocate memory for curve fitting.
+
+	nx = IM_LEN(in, 1)
+
+	call smark (sp)
+	call salloc (x, nx, TY_REAL)
+	call salloc (wts, nx, TY_REAL)
+
+	do i = 1, nx
+	    Memr[x+i-1] = i
+	call amovkr (1., Memr[wts], nx)
+
+	call ic_putr (ic, "xmin", Memr[x])
+	call ic_putr (ic, "xmax", Memr[x+nx-1])
+
+	# If the interactive flag is set then use icg_fit to set the
+	# fitting parameters.  Get_fitline returns EOF when the user
+	# is done.  The weights are reset since the user may delete
+	# points.
+
+	if (interactive) {
+	    call clgstr ("graphics", graphics, SZ_FNAME)
+	    gp = gopen ("stdgraph", NEW_FILE, STDGRAPH)
+
+	    i = strlen (title)
+	    while (lf_getline (ic, gt, in, indata, inline, title)!=EOF) {
+		title[i + 1] = EOS
+	        call icg_fit (ic, gp, "cursor", gt, cv, Memr[x], Memr[indata],
+		    Memr[wts], nx)
+		call amovkr (1., Memr[wts], nx)
+	    }
+	    call gclose (gp)
+	}
+
+	# Loop through each input image line and create an output image line.
+
+	new = YES
+	call amovkl (long(1), inline, IM_MAXDIM)
+	call amovkl (long(1), outline, IM_MAXDIM)
+
+	while (imgnlr (in, indata, inline) != EOF) {
+	    if (impnlr (out, outdata, outline) == EOF)
+		call error (0, "Error writing output image")
+
+	    call ic_fit (ic, cv, Memr[x], Memr[indata], Memr[wts],
+		nx, new, YES, new, new)
+	    new = NO
+
+	    call ic_clean (ic, cv, Memr[x], Memr[indata], Memr[wts], nx)
+
+	    call amovr (Memr[indata], Memr[outdata], nx)
+	}
+
+	call cvfree (cv)
+	call sfree (sp)
+end
+
+
+# LC_IMMAP -- Map images for lineclean.
+
+procedure lc_immap (input, output, in, out)
+
+char	input[ARB]		# Input image
+char	output[ARB]		# Output image
+pointer	in			# Input IMIO pointer
+pointer	out			# Output IMIO pointer
+
+pointer	sp, root, sect
+int	imaccess()
+pointer	immap()
+
+begin
+	# Get the root name and section of the input image.
+
+	call smark (sp)
+	call salloc (root, SZ_FNAME, TY_CHAR)
+	call salloc (sect, SZ_FNAME, TY_CHAR)
+
+	call get_root (input, Memc[root], SZ_FNAME)
+	call get_section (input, Memc[sect], SZ_FNAME)
+
+	# If the output image is not accessible then create it as a new copy
+	# of the full input image.
+
+	if (imaccess (output, 0) == NO)
+	    call img_imcopy (Memc[root], output, false)
+
+	# Map the input and output images.
+
+	in = immap (input, READ_ONLY, 0)
+
+	call sprintf (Memc[root], SZ_FNAME, "%s%s")
+	    call pargstr (output)
+	    call pargstr (Memc[sect])
+	out = immap (Memc[root], READ_WRITE, 0)
+
+	call sfree (sp)
+end
+
+
+# LF_GETLINE -- Get an image line to be fit interactively.  Return EOF
+# when the user enters EOF or CR.  Default line is the first line and
+# the out of bounds lines are silently limited to the nearest in bounds line.
+
+int procedure lf_getline (ic, gt, im, data, v, title)
+
+pointer	ic			# ICFIT pointer
+pointer	gt			# GTOOLS pointer
+pointer	im			# IMIO pointer
+pointer	data			# Image data
+long	v[ARB]			# Image line vector
+char	title[ARB]		# Title
+
+int	i
+char	line[SZ_LINE]
+int	getline(), nscan(), imgnlr()
+
+begin
+	call sprintf (title, SZ_LINE, "%s: Fit line =")
+	    call pargstr (title)
+
+	call printf ("%s ")
+	    call pargstr (title)
+	call flush (STDOUT)
+
+	if (getline(STDIN, line) == EOF)
+	    return (EOF)
+	call sscan (line)
+
+	call amovkl (long (1), v, IM_MAXDIM)
+	do i = 2, max (2, IM_NDIM(im)) {
+	    call gargl (v[i])
+	    if (nscan() == 0)
+		return (EOF)
+	    else if (nscan() != i - 1)
+		break
+
+	    if (IM_NDIM(im) == 1)
+		v[i] = 1
+	    else
+	        v[i] = max (1, min (IM_LEN(im, i), v[i]))
+
+	    call sprintf (title, SZ_LINE, "%s %d")
+		call pargstr (title)
+		call pargl (v[i])
+	}
+
+	call sprintf (title, SZ_LINE, "%s\n%s")
+	    call pargstr (title)
+	    call pargstr (IM_TITLE(im))
+	call ic_pstr (ic, "xlabel", "Column")
+	call gt_sets (gt, GTTITLE, title)
+
+	return (imgnlr (im, data, v))
+end