Initial commit

3 files changed, 1388 insertions, 0 deletions

diff --git a/noao/nproto/ace/doc/detect.hlp b/noao/nproto/ace/doc/detect.hlp
new file mode 100644
index 00000000..ac18c675
--- /dev/null
+++ b/noao/nproto/ace/doc/detect.hlp
@@ -0,0 +1,470 @@
+.help detect Sep00 ace
+.ih
+NAME
+detect -- detect and catalog objects in images
+.ih
+SYNOPSIS
+.ih
+USAGE	
+detect images objmasks catalogs
+.ih
+PARAMETERS
+.ls images
+List of images containing objects to be detected.  The images should generally
+have read and write permission to allow addition of header information.
+However, the task will still run without write access with the consequence
+that the header will not be updated.
+.le
+.ls masks = "!BPM"
+List of bad pixel masks for the images.  This may consist of no bad pixel
+mask specified as the empty string "", a single bad pixel mask to apply to
+all images or a list of bad pixel masks which must match the images list.
+Mask names beginning with "!" are image header keywords which point to the
+bad pixel mask.
+.le
+.ls skys = "SKYFIT"
+List of sky images, constant values, sky fit names, or keyword
+indirection.  If only one value is specified then it applies to all input
+images otherwise the list must match the images list.  Values beginning
+with "!" specify image header keywords containing the image name, constant
+value, or sky fit name to be used.  The value is first checked to see
+if an image with that name exists, then if sky fit keywords are in the
+header, and finally if it is a number.  Sky fit keywords are formed from
+the sky fit name with two digit sequence numbers and are interpreted as
+surface fit coefficients.
+
+If none of these are found then the value is treated as the sky fit name
+to be used to save sky fitting performed by this task.
+.le
+.ls sigmas = "SKYSIG"
+List of sky sigma images, constant values, sigma fit names, or keyword
+indirection.  If only one value is specified then it applies to all input
+images otherwise the list must match the images list.  Values beginning
+with "!" specify image header keywords containing the image name, constant
+value, or sigma fit name to be used.  The value is first checked to see
+if an image with that name exists, then if sigma fit keywords are in the
+header, and finally if it is a number.  Sigma fit keywords are formed from
+the sigma fit name with two digit sequence numbers and are interpreted as
+surface fit coefficients.
+
+If none of these are found then the value is treated as the sigma fit name
+to be used to save sky fitting performed by this task.
+.le
+
+The following parameters specify the output.
+.ls objmasks
+List of output object masks.  If no list is given then no object masks
+will be created.  Otherwise there must be one object mask name for each
+input image.  The object mask name will be recorded in the input image
+header and in any output catalog.
+.le
+.ls catalogs
+List of output catalogs.  If no list is given then no catalogs will be
+created.  Otherwise there must be one catalog name for each input image.
+The catalog name will be recorded in the input image header and in any
+object mask.  The catalog is created as a "table" (see \fBtables\fR
+for information about the tables and general tools to interact with the
+tables).  If the name has an explicit ".fits" extension then a FITS binary
+table is created otherwise an IRAF table (".tab" extension) is created.
+.le
+.ls logfiles = "STDOUT"
+List of output log files.  If no list is given then no output log information
+will be produced.  If only one file is specified it applies to all input
+images otherwise the list of files must match the images list.  Note that
+the special name "STDOUT" corresponds to terminal output.
+.le
+
+The following parameters define the initial sky fit determination.  This is
+only done if no sky image or sky constant value and sigma image or sigma
+constant value are specified.
+# Sky
+.ls newsky = no
+Determine new sky fit if one already exists?  When the specified sky
+corresponds to an existing sky fit (the sky fit coefficients are in the
+image header) then this parameter is used to override that fit with a new
+fit.  Otherwise the fit is used and the initial sky fitting is skipped.
+The sky fitting is also skipped if the specified sky is an image or
+constant.
+.le
+.ls nskylines = 100
+Number of sky sample lines to use.  This number of lines spread evenly
+through the image are used to determine the initial sky fit.
+.le
+.ls skyblk1d = 10
+Sky block size for 1D sky estimation. 
+.le
+.ls skyhclip = 2.
+High sky clipping during 1D sky estimation
+.le
+.ls skylclip = 3.
+Low sky clippling during 1D sky estimation
+.le
+.ls skyxorder = 4
+Sky fitting x order
+.le
+.ls skyyorder = 4
+Sky fitting y order
+.le
+.ls skyxterms = "half" (none|half|full)
+Sky fitting y order
+.le
+
+# Iterated Sky
+.ls skyupdate = no
+Update sky after detection iterations?
+.le
+.ls niterate = 1
+Maximum number of sky iterations
+.le
+.ls skyblk2d = 50
+Sky block size during detection
+.le
+.ls maxskyres = 0.2
+Maximum sky residual for iteration
+.le
+
+# Detection
+.ls convolve = "block 3 3"
+Convolution kernel
+.le
+.ls hsigma = 3.
+Sigma threshold above sky
+.le
+.ls lsigma = 10.
+Sigma threshold below sky
+.le
+.ls hdetect = yes
+Detect objects above sky?
+.le
+.ls ldetect = yes
+Detect objects below sky?
+.le
+.ls minpix = 10
+Minimum number of pixels in detected objects
+.le
+.ls sigavg = 4.
+Sigma of mean flux cutoff
+.le
+.ls sigmax = 4.
+Sigma of maximum pixel
+.le
+.ls bpval = 1
+Output bad pixel value
+.le
+
+# Splitting"
+.ls split = yes
+Split objects?
+.le
+.ls splitmax = INDEF
+Maximum sigma above sky for splitting
+.le
+.ls splitstep = 0.4
+Splitting steps in convolved sigma
+.le
+.ls splitthresh = 5.
+Splitting threshold in sigma
+.le
+.ls sminpix = 10
+Minimum number of pixels in split objects
+.le
+.ls ssigavg = 10.
+Sigma of mean flux cutoff
+.le
+.ls ssigmax = 5.
+Sigma of maximum pixel
+.le
+
+
+# Growing"
+.ls ngrow = 2
+Number of grow rings
+.le
+.ls agrow = 2.
+Area grow factor
+.le
+.ih
+DESCRIPTION
+
+SKY DETERMINATION
+
+A critical part of detecting objects in astronomical images is determining
+the background sky and sky sigma at each point in the image.  In the
+following discussion sky means both the mean sky level and the sky sigma.
+\fBDetect\fR provides for either the user to specify the sky or for the
+task to use a sky fitting algorithm.  The user may specify a sky either as
+another image or as a constant value.  Note that the image name or
+value may be specified either explicitly or with a keyword associated
+with the image.
+
+If the sky is not specified by an image or constant value then a surface
+fit to the sky is used.  The surface fit is recorded in the image header as
+a sequence of keywords with a specified name (the keyword prefix which may
+be up to six characters) and two digit sequence number.  The values of the
+keywords contain the coefficients of the fit.  The the surface fit
+coefficients are defined in the SURFACE FIT section.
+
+Note that it is possible to specify the mean sky and the sky sigma in
+different ways.  When one is given as an image or constant and the other
+as a fit.  The one given as an image or constant will be kept fixed and
+the fit determination and updating will be done only on the other.
+
+The sky surface fit is computed in two stages.  There is an initial
+determination using a subsample of image lines.  Then there is an
+optional update of the sky sample during the object detection step.
+The detection step with sky updating may be iterated a specified number
+of times until the maximum difference in the mean sky is less than some
+amount.
+
+INITIAL SKY DETERMINATION
+
+If an existing surface fit is specified then the parameter \fInewsky\fR
+selects whether a new surface fit is to be computed.  If the value is "no"
+then the initial sky determination is skipped though the detection update
+may still be selected.
+
+The initial sky fit uses a combination of block averaging to reduce the
+number of points in the fitting, one dimensional line fitting with sigma
+clipping rejection to eliminate objects, and finally fitting a two
+dimensional surface to the set of block averages over all the sample lines
+which cover the image.
+
+The parameter \fInskylines\fR defines the number of sample lines across
+the image to be used.  The lines are evenly spaced starting with the
+first line and ending with the last line.  The number of lines affects
+how fast the sky estimation is done.
+
+The pixels from the input line are initially all given unit weight.  Bad
+pixels identified by the input bad pixel mask are excluded by setting their
+weights to zero.  A weighted block average, with the weight of each block
+being the sum of the weights, is computed.  The size of the blocks is given
+by the \fIskyblk1d\fR parameter.  This is done to speed the fitting by
+reducing the number of points.  Note that when all pixels in a block have
+zero weight due to the bad pixel mask or subsequent rejection the weight of
+the composite block average point is zero.
+
+If only one of sky mean and sky sigma quantities is being determined with
+the other quantity given by an input image, constant, or previous fit
+then those values are simple block averaged with the same block size
+to produce sample points for the mean sky or sky sigma.  Note that the
+sky sigma of the sample points also requires division by the square root
+of the block size to give the sky sigma per block average point.  The
+line fitting described next is then skipped for this quantity.
+
+The weighted one dimensional line fitting to the block averages uses
+Chebyshev polynomials of order given by the \fIskyxorder\fR.  Note that
+this order is the number of polynomial terms, which is one higher than the
+maximum power of the polynomial so that a value of 3 corresponds to a
+quadratic polynomial.
+
+When the mean sky is being determined, the line fitting is performed and
+the fitted values at the block centers are evaluated.
+
+When the sky sigma is being determined, the absolute value of the residuals
+relative to the mean sky divided by 0.7979 are computed.  A gaussian noise
+distribution will have a mean value of this quantity equal to the sigma of
+the distribution.  In other words, the mean of the absolute deviations of a
+gaussian distribution is 0.7979 times sigma.  By fitting a function to
+these residual values a position variable estimate of the sky sigma is
+obtained without needing to compute standard deviations over some set of
+points.  The fitted values at the block centers are evaluated to give the
+sky sigmas for the block averaged data.
+
+With the set of block averaged data points and estimated mean skys and sky
+sigmas points that deviate by more than the number of sigma given by the
+\fIskyhclip\fR and \fIskylclip\fR parameters are rejected by setting their
+weights to zero.  The line fitting is then repeated until no points are
+rejected with a maximum of 10 iterations.
+
+When the iteration completes the block average points for that image line
+are accumulated for a two dimensional surface fit.  Note that the weights
+are used to exclude rejected averages and to weight blocks that had fewer
+points due to bad pixels.  The surface fit is a two dimensional Chebyshev
+polynomial of orders given by the \fIskyxorder\fR and \fIskyyorder\fR.  The
+orders have the same meaning as in the one dimensional polynomial, namely
+the number of terms in powers of x and y.  There are also cross terms which
+are a mixture of powers of both x and y.  The \fIskyxterms\fR select
+whether to use any cross terms, only cross terms whose total power does not
+exceed the maximum of the pure x and y terms, or all combinations of
+powers.
+
+After all the sample lines are completed the final surface fits are
+computed.  The coefficients of the fits are written to the image header
+under the specified sky fit names and the fits are passed on to the
+detection phase.  Note that if the input image is read only then the
+fit will not be written to the header but the task continues.
+
+UPDATED TO SKY DURING DETECTION
+
+
+DETECTION
+
+The detection of objects in an image is conceptually quite simple.  Each
+pixel is compared against the expected sky at that point and if it is
+more that a specified number of sky sigma above the sky it is a candidate
+object pixels.  Candidate object pixels are grouped into objects on the basis
+of being connected along the eight neighboring directions.  The candidate
+object is then accepted if it satisfies the criteria of a minimum
+number of pixels, a sufficiently significant maximum pixel, and a sufficiently
+significant flux above sky.
+
+To detect faint objects where individual pixels are not significantly above
+the sky but all pixels taken together are significant a detection filter is
+applied.  This consists of applying a convolution function to the image and
+performing the detection described in the previous paragraph on the
+convolved pixels with the sky sigma suitable adjusted for the convolution.
+The convolution acts as an optimizing filter for objects with shapes
+corresponding to the convolution weights.  The remaining discussion
+is in terms of the convolved pixel values.  The case of no convolution
+can be thought of as a convolution with a delta function though the
+implementation is not done as a convolution for efficiency.
+
+Two other options to the detection are to also find pixels that are
+significantly below sky (using an independent threshold to that used for
+detecting pixels above sky) and form them into "dark" objects and to
+take the remaining pixels that are not significantly above or below the
+sky and use them to define a sky sample for output or for updating the
+initial sky.
+
+We now go into more detail.  The background sky and sky sigma against which
+the detection is performed is initially set as described earlier.  If desired
+the sky pixels may be accumulated to update the sky.  After updating the
+sky the detection step may be repeated using the new sky.  This is
+discussed  futher when we reach the end of the detection step description.
+
+The convolution is specified by the \fIconvolve\fR parameter.  The values for
+this parameter and the definition of the convolution are given in the
+CONVOLUTION DETECTION FILTER section.  The input pixel data is convolved
+and the sky sigma is appropriately adjusted.
+
+When the central pixel in the convolution is flagged as a bad pixel by the
+bad pixel mask (any non-zero value is a bad pixels) then the convolved
+value is considered to be a bad pixel.  If an output object masks is
+specified the pixel will be marked  with the value specified by the
+\fIbpval\fR parameter.  The value may be set to not show the bad pixel in
+the object mask, to set all input bad pixels to some value, or to pass the
+input bad pixel value to the object mask.  Note that bad pixel masks in the
+object mask must be between 1 and 10 to avoid confusion with the values
+used to identify objects.  If other pixels in the convolution are flagged
+as bad pixels they are excluded from the convolution and the
+convolved sky sigma is adjusted but the convolution value is still used
+as a valid image pixel for detection.
+
+The sigma threshold for pixels to be detected as part of an object above
+sky is given by the \fIhsigma\fR.  This number is multiplied by the sky
+sigma to get the deviation from sky.  As noted earlier the sky sigma is
+for the convolved pixels and the 
+
+
+CONVOLUTION DETECTION FILTER
+
+
+The convolution detection filter is specified with the \fIconvolve\fR
+parameter.  There is only one convolution that can be specified and it applies
+to all input images in a list.  If a null string ("") is specified
+then no convolution is performed.  The task has been optimizations for
+this case to avoid treating this as a 1x1 convolution and to avoid extra
+memory allocations required when a convolution is done.
+
+The convolved value at pixel (i,j), denoted I(i,j), within an image of size
+CxL is defined by
+
+.nf
+    I_convolved(i,j) = sum_kl{I_unconvolved(m,n)*W(k,l)} / sum_kl{W(k,l)}
+.fi
+
+where I(m,n) is the unconvolved value at pixel (m,n), W(k,l) are the NX x
+NY (both must be odd) convolution weights, sum_kl is the double sum over k
+and l, and
+
+.nf
+    m' = i + k - (NX+1)/2	for k = 1 to NX
+    n' = j + l - (NY+1)/2	for l = 1 to NY
+
+    m = m' (1<=m'<=C)	m = 1-m' (m'<1)	  m = 2C-m' (m'>C)
+    n = n' (1<=n'<=L)	n = 1-n' (n'<1)	  n = 2L-n' (m'>L)
+.fi
+
+The last two lines represent boundary reflection at the edges of the image.
+
+The sky sigma of a convolved pixel is approximated by
+
+.nf
+    sigma_convolved(i,j) = sigma_unconvolved(i,j) / sum_kl{W(k,l)}
+.fi
+
+In the presence of bad pixels identified by a bad pixel mask the convolution
+weight applied to a bad pixel is set to zero.  The sum of the weights
+used to normalize the convolution is then modified from the situation with
+no bad pixels.  This will correct the convolved pixel value for the missing
+data and the estimated sky sigma is appropriately larger.
+
+A convolution can be computational slow, especially for larger sizes.
+The implementation of the convolution has been optimized to recognize
+bilinear symmetries or lines which are scaled versions of other lines.
+So if possible such symmetries should be used.  The "block", "bilinear",
+and "gauss" special convolutions described below have such symmetries.
+
+There is also an overhead in checking for bad pixels.  The convolution
+has an optimization to avoid such checks in the case where no bad pixel
+mask is specified.
+
+The \fIconvolve\fR parameter is a string which can take one of the
+following forms.
+
+.ls ""
+There is no convolution or, equivalently, NX=1, NY=1.
+.le
+.ls @[filename]
+The weights are given in the specified file.  The format consists of lines
+of whitespace separated values.  The number of values on each line must be
+the same and defines NX and the number of lines defines NY.
+.le
+.ls block [NX] [NY]
+The weights are all the same and the convolution size is given by the
+two numbers following the word "block".
+.le
+.ls bilinear [NX] [NY]
+The weights are the bilinear matrix product of triangular one dimensional
+matrices of sizes given by the two numbers following the word "bilinear".
+The weights are described by the matrix product relation 
+
+.nf
+    [1 ... (NX+1)/2 ... 1] * Transpose{[1 ... (NY+2)/2 ... 1]}
+.fi
+
+For example for NX=5, and NY=3 the weights would be
+
+.nf
+    1 2 3 2 1
+    2 4 6 4 2
+    1 2 3 2 1
+.fi
+.le
+.ls gauss [NX] [NY] [SX] [SY]
+The weights are bidimensional gaussian values on a grid of size NX by NY
+with sigma values SX and SY (real numbers) in units of pixel spacing.
+.le
+.ls [W(1,1)] ... [W(NX,1)], ..., [W(1,NY)] ... [W(NX,NY)]
+The weights are specified as a string of real values.  The values are
+whitespace separated within each line and the lines are delimited by
+comma.  For example
+
+.nf
+                               1 2 1
+    1 2 1, 2 3 2, 1 2 1  ==>   2 3 2
+                               1 2 1
+.le
+
+When a logfile is defined the weights are included in the log output.
+
+
+OBJECT MASKS
+
+.ih
+EXAMPLES
+.ih
+REVISIONS
+.ih
+SEE ALSO
+.endhelp
diff --git a/noao/nproto/ace/doc/installation.hlp b/noao/nproto/ace/doc/installation.hlp
new file mode 100644
index 00000000..c399ad4c
--- /dev/null
+++ b/noao/nproto/ace/doc/installation.hlp
@@ -0,0 +1,208 @@
+.help installation Jan01 ace
+
+.ce
+\fBACE: Astronomical Cataloging Environment\fR
+.ce
+Release Notes and Installation Instructions
+
+.sh
+SUMMARY
+The ACE external package is used to catalog objects in images and manipulate
+the catalogs.
+
+.sh
+RELEASE INFORMATION
+The following summary only highlights the major changes.  There will also
+be minor changes and bug fixes.
+
+.ls V0.2: January 27, 2001
+Alpha test version.
+.le
+.sh
+INSTALLATION INSTRUCTIONS
+Installation of this external package consists of obtaining the files,
+creating a directory containing the package, compiling the executables or
+installing precompiled executables, and defining the environment to load
+and run the package.  The package may be
+installed for a site or as a personal installation.  If you need help with
+these installation instructions contact iraf@noao.edu or call the IRAF
+HOTLINE at 520-318-8160.
+.ls [arch]
+In the following steps you will need to know the IRAF architecture
+identifier for your IRAF installation.  This identifier is similar to the
+host operating system type.  The identifiers are things like "ssun" for
+Solaris, "alpha" for Dec Alpha, and "linux" or "redhat" for most Linux
+systems.  The IRAF architecture identifier is defined when you run IRAF.
+Start the CL and then type
+
+.nf
+    cl> show arch
+    .ssun
+.fi
+
+This is the value you need to know without the leading '.'; i.e. the
+IRAF architecture is "ssun" in the above example.
+.le
+.ls [1-site]
+If you are installing the package for site use, login as IRAF
+and edit the IRAF file defining the packages.
+
+.nf
+    % cd $hlib
+.fi
+
+Define the environment variable ace to be the pathname to
+the ace package root directory.  The '$'
+character must be escaped in the VMS pathname and UNIX pathnames must be
+terminated with a '/'.  Edit extern.pkg to include the following.
+
+.nf
+    reset ace = /local/ace/
+    task  ace.pkg = ace$ace.cl
+.fi
+
+Near the end of the hlib$extern.pkg file, update the definition of
+helpdb so it includes the ace help database, copying the syntax
+already used in the string.  Add this line before the line 
+containing a closing quote:
+
+.nf
+    ,ace$lib/helpdb.mip\
+.fi
+.le
+.ls [1-personal]
+If you are installing the package for personal use define a host
+environment variable with the pathname of the directory where the package
+will be located (needed in order to build the package from the source
+code).  Note that pathnames must end with '/'.  For example:
+
+.nf
+    % setenv ace /local/ace/
+.fi
+
+In your login.cl or loginuser.cl file make the following definitions
+somewhere before the "keep" statement.
+
+.nf
+    reset ace = /local/ace/
+    task  ace.pkg = ace$ace.cl
+    printf ("reset helpdb=%s,ace$lib/helpdb.mip\nkeep\n",
+        envget("helpdb")) | cl
+    flpr
+.fi
+
+If you will be compiling the package, as opposed to installing a binary
+distribution, then you need to define various environment variables.
+The following is for Unix/csh which is the main supported environment.
+
+.nf
+    # Example
+    % setenv iraf /iraf/iraf/             # Path to IRAF root (example)
+    % source $iraf/unix/hlib/irafuser.csh # Define rest of environment
+    % setenv IRAFARCH ssun                # IRAF architecture
+.fi
+
+where you need to supply the appropriate path to the IRAF installation root
+in the first step and the IRAF architecture identifier for your machine
+in the last step.
+.le
+.ls [2]
+Login into IRAF.  Create a directory to contain the package files and the
+instrument database files.  These directory should be outside the standard
+IRAF directory tree.
+
+.nf
+    cl> mkdir ace$
+    cl> cd ace
+.fi
+.le
+.ls [3]
+The package is distributed as a tar archive for the
+sources and, as an optional convenience, a tar archive of the executables
+for select host computers.  Note that IRAF includes a tar reader.  The tar
+file(s) are most commonly obtained via anonymous ftp.  Below is an example
+from a Unix machine where the compressed files have the ".Z" extension.
+Files with ".gz" or ".tgz" can be handled similarly.
+
+.nf
+    cl> ftp iraf.noao.edu (140.252.1.1)
+    login: anonymous
+    password: [your email address]
+    ftp> cd iraf/extern
+    ftp> get ace.readme
+    ftp> binary
+    ftp> get ace.tar.Z
+    ftp> get ace-bin.<arch>.Z  (optional)
+    ftp> quit
+    cl> !uncompress ace.tar
+    cl> !uncompress ace-bin.<arch> (optional)
+.fi
+
+The readme file contains these instructions.  The <arch> in the
+optional executable distribution is replaced by the IRAF architecture
+identification for your computer.
+
+Upon request the tar file(s) may be otained on tape for a service
+charge.  In this case you would mount the tape use rtar to extract
+the tar files.
+.le
+.ls [4]
+Extract the source files from the tar archive using 'rtar".
+
+.nf
+    cl> softools
+    so> rtar -xrf ace.tar
+    so> bye
+.fi
+
+On some systems, an error message will appear ("Copy 'bin.generic'
+to './bin fails") which can be ignored.
+Sites should leave the symbolic link 'bin' in the package root
+directory pointing to 'bin.generic' but can delete any of the
+bin.<arch> directories that won't be used.  If there is no binary
+directory for the system you are installing it will be created
+when the package is compiled later or when the binaries are installed.
+
+If the binary executables have been obtained these are now extracted
+into the appropriate bin.<arch> directory.
+
+.nf
+    # Example of sparc installation.
+    cl> cd ace
+    cl> rtar -xrf ace-bin.sparc      # Creates bin.sparc directory
+.fi
+
+The various tar file can be deleted once they have been
+successfully installed.
+.ls [5]
+For a source installation you now have to build the package
+executable(s).  The "tables" package must be installed first if not
+already available.  First you configure the package for the particular
+architecture.
+
+.nf
+    cl> cd ace
+    cl> mkpkg <arch>            # Substitute sparc, ssun, alpha, etc.
+.fi
+
+This will change the bin link from bin.generic to bin.<arch>.  The binary
+directory will be created if not present.  If an error occurs in setting
+the architecture then you may need to add an entry to the file "mkpkg".
+Just follow the examples in the file.
+
+To create the executables and move them to the binary directory
+
+.nf
+    cl> mkpkg -p ace 		# build executables
+    cl> mkpkg generic           # optionally restore generic setting
+.fi
+
+Check for errors.  If the executables are not moved to the binary directory
+then step [1] to define the path for the package was not done correctly.
+The last step restores the package to a generic configuration.  This is not
+necessary if you will only have one architecture for the package.
+.le
+
+This should complete the installation.  You can now load the package
+and begin testing and use.
+.endhelp
diff --git a/noao/nproto/ace/doc/objmasks.hlp b/noao/nproto/ace/doc/objmasks.hlp
new file mode 100644
index 00000000..1c35c4c9
--- /dev/null
+++ b/noao/nproto/ace/doc/objmasks.hlp
@@ -0,0 +1,710 @@
+.help objmasks Jan02 nproto
+.ih
+NAME
+objmasks -- detect objects in images and create masks and sky maps
+.ih
+SYNOPSIS
+.ih
+USAGE	
+objmasks images objmasks skys
+.ih
+PARAMETERS
+.ls images
+List of images or multiextension files for which object masks are desired.
+.le
+.ls objmasks
+List of object masks to be created.  This list must match the input list.
+Multiextension input files will produce multiextension mask files.  If the
+input image is writable, the name of the created mask will recorded in the
+image header.  Note that it is possible to specify a null image to
+not produce an output mask.  This might be done if the background sky
+or sky sigma maps are desired or to just see the log information.
+.le
+
+.ls omtype = "numbers" (boolean|numbers|colors|all)
+The type of encoding for the object mask values.  In all cases non-object pixels
+(that is background) have mask values of zero.  The choices for the mask
+values are "boolean", "numbers", "colors", and "all".  These are described
+in the \fIOutput Data\fR section.
+.le
+.ls skys = "", sigmas = ""
+Optional lists of input or output sky and sigma maps.  Maps are either
+constant values or images which are interpolated to the size of the input
+images.  If a list is given it must match the input \fIimages\fR list.
+If constant values or existing maps are specified then those are used
+without change.  If a new filename is given then an output file is created
+with the values computed by the task.  Multiextension input images create
+or apply the same extension names to the specified sky or sigma files.
+Constant input values apply to all extensions.  The sigma values are
+per single input image pixel.
+.le
+.ls masks = "!BPM"
+List of bad pixel masks for the input images.  Non-zero masks values are
+ignored in the object detection and are passed on to the output object
+masks based on the \fIomtype\fR parameter.  An empty list applies no bad
+pixel mask, a single mask applies to all input images, and a matching
+list matches the masks with the input image.  A mask is specified by a
+filename or by reference to a filename given by the value of a header
+keyword in the input image.  A header keyword reference is made with the
+syntax "!<keyword>" where <keyword> is the desired keyword with case
+ignored.  For multiextension files the input masks may be either a
+multiextension file with matching extension names or a directory of
+pixel list files with the extension names as filenames.
+.le
+.ls extnames = ""
+Extensions to select from multiextension files.  A null string matches all
+extension names.  Otherwise the parameter is a comma separated list of
+patterns that match the entire extension name.  Thus, an explicit list of
+extension names may be specified or the pattern matching characters '?' for
+any character or '[]' for a set of characters may be used.  The set may
+include ranges in ascii order by using hyphens; i.e. 1-3 matches the
+characters 1, 2, and 3.
+.le
+.ls logfiles = "STDOUT"
+List of output log files.  If no list is given then no output log information
+will be produced.  If only one file is specified it applies to all input
+images otherwise the list of files must match the images list.  Note that
+the special name "STDOUT" corresponds to terminal output.
+.le
+
+.ls blkstep = 1
+The mean and sigma of the background or sky pixels are determined in a
+first pass through the image.  If \fIblkstep\fR is one all lines are used.
+To skip lines in order to speed up this computation, the parameter may be
+set to a larger value to define the increment between lines.  However, the
+task will enforce a preset minimum number to insure a sufficient sample.
+.le
+.ls blksize = -10
+The background mean sky and sky sigma are determined in a set of square
+blocks from which the values are linearly interpolated to each point in the
+input image.  The size of the blocks may be specified as a number of blocks
+spanning the smaller image dimension by using a negative integer value.
+Or the size may be specified as the number of pixels across a block.
+The task will enforce a preset minimum number of pixels per block which may
+require using bigger blocks than specified.  The background determination
+algorithm is described further in the "Background Determination" section.
+.le
+
+.ls convolve = "block 3 3"
+Convolution filter to be applied prior to threshold detection.  The
+convolution filter is defined by a set of weights in a 2D array.  These
+may be specified in files or with certain forms given by special strings.
+The options are described in the "Convolution Filter" section.
+.le
+.ls hsigma = 3., lsigma = 10.
+Object pixels are identified by sigma thresholds about the mean background
+based on the estimated background sigma at each point in the image.
+The sigma factors are specified in terms of the "per pixel" sigma before
+convolution.  The \fIhsigma\fR value is the "high" or above background
+limit and the \fIlsigma\fR value is the "low" or below background limit.
+Typically detections are one-sided, such as detecting objects above
+the background, and so the thresholds need not be equal.
+.le
+.ls hdetect = yes, ldetect = no
+Identify objects as pixels which are above the background (\fIhdetect\fR)
+and below the background (\fIldetect\fR)?  If objects are detected but the
+corresponding parameter is no then the output mask will not include those
+objects.
+.le
+.ls neighbors = "8" (8|4)
+The threshold selected pixels are associated with other neighboring pixels to
+form an object.  The criterion for a neighbor being part of the
+same object is defined by this parameter.  The choices are "8" for
+pixels touching in any of the 8 directions or "4" to identify neighbors
+as only horizontal or vertically adjacent.
+.le
+.ls minpix = 6
+The minimum number of neighboring pixels which define an acceptable object.
+.le
+.ls ngrow = 2, agrow = 2.
+After an object is identified as a set of threshold detected pixels,
+additional neighboring pixels may be added to the object.  This allows
+expanding the object into the faint wings of the light distribution.  The
+additional pixels are those which touch the boundary pixels.  Pixels are
+added in multiple passes, each time extending the previous boundary.  The
+parameter \fIngrow\fR (an integer value) defines the maximum number of
+boundary extensions.  The parameter \fIagrow\fR (a real value) specifies
+the maximum increase in area (number of pixels) from the original
+detection.
+.le
+.ih
+DESCRIPTION
+\fBOBJMASKS\fR is a task for creating masks covering objects in images.
+An optional secondary product of this task is to produce background
+and sigma maps.  Objects are identified by threshold sigma detection.
+These object masks may be used by other applications to exclude the object
+data or focus on the objects.  The detection consists of determining a
+smooth, spatially variable mean background and background sigma (if no
+input maps are provided), convolving the data by an optional filter to
+optimize detection of faint sources, collecting pixels satisfying the
+detection thresholds, assigning neighboring pixels to a common object,
+applying a minimum number of pixels test to the objects, and growing
+objects to extend into the wings of the object light distribution.
+The last step is writing out the identified object pixels as a mask.
+
+1. Input Data
+
+The input data consists of one or more 2D images.  The images are assumed
+to  contain a moderately smooth background and multiple sources or
+objects.  This task is most useful for images with large numbers of small
+sources rather than one large object such as a nearby galaxy.  The input
+images, specified by the \fIimages\fR parameter, may be individual images
+(which includes images selected from multiextension files as explicit
+image extensions) or multiextension files specified by a root filename.  In
+the latter case the image extension names selected by the \fIextnames\fR
+parameter are used.
+
+Background means and sigmas (specified per image pixels) may be specified
+by "maps".  These may be constant numerical values or images.  The map
+images will be linearly interpolated to the size of the input images.
+For multi-extension input data, constant map values apply to all extensions
+and maps are also multiextension files with map images having the same
+extension names.
+
+Bad pixel masks may be associated with the input images to
+exclude pixels from the background and object determinations.  These
+bad pixels are also included in the output object masks.  The bad pixel
+masks are specified by the \fImasks\fR parameter.  This parameter may
+identify a mask by a filename or a keyword.  A single mask may be
+specified to apply to all images or a matching list of masks may be
+given.
+
+The masks are in one of the supported mask formats.  As of IRAF V2.12 this
+includes pixel list (.pl) files and FITS "type=mask" extensions.  When the
+input files are multiextension files, the selected extension names are
+appended to the specified mask filename to select masks with the same
+extension name.  If a mask file of the form "name[ext]" is not found
+the task will treat the filename as a directory of pixel list files and
+select the pixel list file with the extension name; i.e. "name/ext.pl".
+
+2. Output Data
+
+The output of this task are object masks, sky maps, sigma maps, and log
+information.  The output object masks default to mask type extensions.  If an
+extension name is not specified explicitly the default extension name
+"pl" is created.  To select a pixel list output format an explicit ".pl"
+extension must be used.
+
+When the input data are multiextension files, the output masks, mean sky
+maps, and sky sigma maps will be multiextension files with the specified
+rootnames and the same extension name as the input.
+
+The output mask values identify non-object pixels with zero.  The non-zero
+values are encoded as selected by the \fIomtype\fR parameter.  The choices
+are:
+
+.ls "boolean"
+All object and bad pixels have a mask value of one; i.e. the output masks
+consists only of the values 0 and 1.
+.le
+.ls "numbers"
+Input bad pixels values between 1 and 10 preserve their value and all
+other input mask values are mapped to 10.  The object mask pixels have
+object numbers starting with 11.  The object numbers are assigned by
+the task (roughly in order from the first line to the last line) and
+all pixels from a single object have the same unique object number.
+.le
+.ls "colors"
+Input bad pixels are mapped to output values of one.  The object numbers
+are modulo 8 plus 2; i.e. values between 2 and 9.  The purpose of this
+numbering is to allow mapping to the nine standard display colors for an
+interesting overlay with the \fBdisplay\fR task and "ocolors='+203'".
+.le
+.ls "all"
+This is the same as "numbers" except that bits 24 to 27 in the mask values
+are used for various purposes.  In particular bit 24 is set for the boundary
+pixels.  This numbering will be used in the future by special tasks.
+.le
+
+Output mean sky and sky sigma maps consist of the mean and sigma values
+in blocks as described in the "Background Determination" section.
+Therefore, the size of the map images are smaller than the input data images.
+These maps need to be interpolated to the size of the input image
+to obtain the values used for particular pixels in the data images.
+This interpolation expansion is done automatically by some tasks such
+as \fBmscred.rmfringe\fR.
+
+The log output provides information about the files, the phase of the
+processing, some of the parameters, and the convolution filter weights.
+The output begins with the task identifier ACE.  This is because this
+prototype task is a first release piece of a major package called ACE
+(Astronomical Cataloging Environment), which is under development.
+
+3. Background Determination
+
+Detection of sources in an image begins with determining the background.
+By this we mean estimating the probability distribution of the background
+pixel values at every pixel in the image.  In practice we only estimate
+the central value and width and assume a normal distribution for evaluating
+the significance of deviations from the central value.  Since we normally
+won't have a sample of values at each pixel the distribution is
+determined from a sample of nearby pixels.
+
+In this discussion the central value of a distribution is denoted by <I>.
+It is estimated by the mean or mode of the sample.  The width of the
+distribution about <I> is denoted by <S> and is estimated by the absolute
+mean residual converted to the standard deviation of a normal distribution
+with the same absolute mean residual.  The normal deviation of a value I
+from the distribution is defined as R = (I - <I>) / <S>.
+
+The background may be specified by input maps for one or both of the
+background quantities.  The maps may be constant values which apply
+to all pixels or a grid of values given in an image which are linearly
+interpolated to the full size of the input data.  For those quantities
+which are not input the following algorithm is used for computing
+a map.  The maps may be output and used as a product of this task.
+
+The background and/or sigma are estimated in two initial passes through the
+data.  The first pass algorithm fits linear functions to a subsample of
+lines using sigma clipping iteration to eliminate objects.  The subsample
+is used to speed up the algorithm and is reasonable since only linear
+functions are used.  Each sample line is block averaged in blocks of 10
+pixels and a linear function is fit by least squares to obtain an estimate
+for <I> along the line.  The fitting weights are the number of good pixels
+in each block average after elimination of bad pixels specified by the
+user in a bad pixel mask.  The absolute values of the residuals are also
+fit to produce a constant function for <S>.
+
+To exclude objects from affecting these estimates the fitting is iterated
+using sigma clipping rejection on the normal deviations R.  In the
+first iteration the fitting function for <S> is a constant and in
+subsequent steps a linear fit is used.  When the sigma clipping iteration
+rejects no more data, the remaining block averages, absolute residuals, and
+weights are used to fit a 2D plane for both <I> and <S>.  The <S> surface
+is a constant in order to avoid potential negative sigma values.
+
+This first pass algorithm is fast and produces good estimates for the
+planar approximation to the background.  The second pass divides the image
+into large, equal sized blocks, as specified by the \fIblksize\fR
+parameter, and estimates <I> and <S> in each block.  The size of the blocks
+needs to be large enough to give good estimates of the statistics though
+small enough to handle the scale of variations in the sky.  Each block is
+divided into four subblocks for independent estimates which are then
+combined into a final value for the block.  As with the first pass, the
+second pass can be speeded up by using a subsample of lines (parameter
+\fBblkstep\fR) provided some minimum number of lines per subblock is
+maintained.
+
+The background estimates in each subblock are made using histograms of the
+normal deviations R computed relative to the first pass estimates of <I>
+and <S>.  When pixels are added into the histogram the <I> and <S> used to
+compute R are accumulated into means of these quantities in order
+to convert estimates from the normalized deviation histogram back into data
+values.  The histograms are truncated at +/-2.5 and have bin widths
+determined by requiring a specified average bin population based on the
+number of pixels in the block.  Typically the bin population is of order
+500.  The histogram truncation is essentially an object-background
+discrimination.
+
+When all the pixels in a subblock have been accumulated, new estimates of
+<I> and <S> are computed.  If the number of pixels in the histogram is
+less than two-thirds of the subblock pixels the estimates are set to be
+indefinite.  This flags the subblock as too contaminated by objects to be
+used.  All subblock neighbors, which may cross the full block boundaries,
+are also rejected to minimize contamination by the wings of big galaxies
+and very bright stars.
+
+If the histogram has enough pixels, the bin populations are squared to
+emphasize the peak of the distribution and reduce the effects of the
+truncated edges of the histogram.  Because of noise and the fine binning of
+the histogram, a simple mode cannot be used and squaring the bin numbers
+helps to approach the mode with a centroid.  Squaring the bin values and
+then computing the centroid can also be thought of as a weighted centroid.
+
+Generally a mode is considered the best estimate to use for the central
+value <I> of the sky distribution.  But it is unclear how to best estimate
+the mode without an infinite number of pixels.  One could do something like
+fit a parabola to the histogram peak.  But instead we use the empirical
+relation for a skewed distribution between the mean, mode, and median;
+<I>=mean-3*(mean-median).  The mean is the weighted centroid and the median
+is obtained numerically from the histogram using linear interpolation to
+get a subbin value.
+
+The <S> values are obtained from the absolute mean residual of the
+unweighted histogram about the previously derived central value <I> of the
+histogram.  The conversion to a standard deviation is made by computing the
+ratio between the standard deviation and mean absolute deviation of a
+Gaussian distribution.  The standard value over the entire distribution
+cannot be used because the histogram is truncated.  However, it is easy to
+numerically compute the ratio with the same truncation.
+
+Once <I> and <S> are obtained in bin numbers it is converted to data
+values by using the mean and sigma of the input pixel values used
+to create the histogram.
+
+The averages of the subblock <I> and <S> values which are not indeterminate
+in each block are computed.  If any of the full blocks are indeterminate
+when all the subblocks have been eliminated as contaminated, values are
+obtained for them by interpolation from nearby blocks.  The block values
+are then linearly interpolated to get background values for every
+pixel in the input image.
+
+Note that the background pixels used in the block algorithm before
+detection are derived by simple sigma clipping of the histogram values
+around the planar background.  If an output map for either the mean
+values or the sigmas is specified then during the object detection stage
+the background and sigmas are updated using the detected sky pixels about
+the initial block sampled background.  This is a more sensitive selection
+of sky pixels since convolution filtering can exclude pixels from faint
+objects and the wings of all objects.  The new set of sky pixels are
+accumulated and used in the same way as described earlier.
+
+4. Convolution Filters
+
+In order to improve the detection of faint sources dominated by the
+background noise, the input data may be convolved to produce filtered
+values in which the noise has been suppressed.  The threshold detection
+is then performed on the filtered data values.
+
+The convolution detection filter is specified with the \fIconvolve\fR
+parameter.  There is only one convolution that can be specified and it
+applies to all input images in a list.  If a null string ("") is specified
+then no convolution is performed.  The task has been optimizations for this
+case to avoid treating this as a 1x1 convolution and to avoid extra memory
+allocations required when a convolution is done.
+
+The convolved value at pixel (i,j), denoted I'(i,j), is defined by
+
+.nf
+    I'(i,j) = sum_kl{I(m,n)*W(k,l)} / sum_kl{W(k,l)}
+.fi
+
+where I(m,n) is the unconvolved value at pixel (m,n), W(k,l) are the NX x
+NY (both must be odd) convolution weights, sum_kl is the double sum over k
+and l, and
+
+.nf
+    m' = i + k - (NX+1)/2	for k = 1 to NX
+    n' = j + l - (NY+1)/2	for l = 1 to NY
+
+    m = m' (1<=m'<=C)	m = 1-m' (m'<1)	  m = 2C-m' (m'>C)
+    n = n' (1<=n'<=L)	n = 1-n' (n'<1)	  n = 2L-n' (m'>L)
+.fi
+
+The size of the image is C x L.  The last two lines represent boundary
+reflection at the edges of the image.
+
+The sky sigma of a convolved pixel is approximated by
+
+.nf
+    sigma'(i,j) = sigma(i,j) / sum_kl{W(k,l)}
+.fi
+
+In the presence of bad pixels specified in the bad pixel mask the
+convolution weight applied to a bad pixel is set to zero.  If the central
+pixel is bad then the convolved value is also considered to be bad.  The
+sum of the weights used to normalize the convolution is then modified from
+the situation with no bad pixels.  This will correct the convolved pixel
+value for the missing data and the estimated sky sigma is appropriately
+larger.  Since there is an overhead in checking for bad pixels the
+convolution has an optimization to avoid such checks in the case where no
+bad pixel mask is specified.
+
+A convolution can be computational slow, especially for larger convolution
+kernel sizes.  The implementation of the convolution has been optimized to
+recognize bilinear symmetries or lines which are scaled versions of other
+lines.  So if possible users should chose convolutions with such symmetries
+to be most efficient.  The "block", "bilinear", and "gauss" special
+convolutions described below all have such symmetries.
+
+The \fIconvolve\fR parameter is a string with one of the following forms.
+
+.ls ""    
+There is no convolution or, equivalently, NX=1, NY=1.
+.le
+.ls @[filename]
+The weights are given in the specified file.  The format consists of lines
+of whitespace separated values.  The number of values on each line must be
+the same and defines NX and the number of lines defines NY.
+.le
+.ls block [NX] [NY]
+The weights are all the same and the convolution size is given by the
+two numbers following the word "block".  This is a moving block average
+filter.
+.le
+.ls bilinear [NX] [NY]
+The weights are the bilinear matrix product of triangular one dimensional
+matrices of sizes given by the two numbers following the word "bilinear".
+The weights are described by the matrix product relation 
+
+.nf
+    [1 ... (NX+1)/2 ... 1] * Transpose{[1 ... (NY+2)/2 ... 1]}
+.fi
+
+For example for NX=5, and NY=3 the weights would be
+
+.nf
+    1 2 3 2 1
+    2 4 6 4 2
+    1 2 3 2 1
+.fi
+.le
+.ls gauss [NX] [NY] [SX] [SY]
+The weights are bidimensional gaussian values on a grid of size NX by NY
+with sigma values SX and SY (real numbers) in units of pixel spacing.
+.le
+.ls [W(1,1)] ... [W(NX,1)], ..., [W(1,NY)] ... [W(NX,NY)]
+The weights are specified as a string of real values.  The values are
+whitespace separated within each line and the lines are delimited by
+comma.  For example
+
+.nf
+                               1 2 1
+    1 2 1, 2 3 2, 1 2 1  ==>   2 3 2
+                               1 2 1
+.fi
+.le
+
+When a logfile is defined the convolution weights are included in the
+output.
+
+5. Object Detection
+
+The detection of objects in an image is conceptually quite simple once the
+background is known.  If an input pixel, before any convolution, is
+identified in the bad pixel mask the output object mask pixel is also
+identified as bad.  Otherwise the input data is convolved as described
+previously.
+
+Each convolved pixel is compared against the expected background at that
+point and, if it is more that a specified number of convolution adjusted
+background sigma above (\fIhsigma\fR) or below (\fIlsigma\fR) the
+background, it is identified as a candidate object pixel.  Candidate object
+pixels, with the same sense of deviation, are grouped into objects on
+the basis of being connected along the four or eight neighboring directions
+as specified by the \fIneighbor\fR parameter.  The candidate object is then
+accepted if it satisfies the minimum number of pixels (\fIminpix\fR) in
+an object and the \fIhdetect\fR or \fIldetect\fR parameter selects that
+type of object.  The accepted objects are assigned sequential numbers
+beginning with 11.  The object numbers are used, as described in the
+section on the output data, to set the output object mask values.
+
+If an output mean sky or sigma map is requested, the output is that
+updated by the sky pixels identified during the detection.
+
+6. Object Growing
+
+Astronomical objects do not have sharp edges but have light distributions
+that merge into the background.  This is due not only to the nature of
+extended sources but to the atmospheric and instrument point spread function
+effects on unresolved sources.  In order to include pixels which extend
+away from the threshold detection and contain some amount of light
+apart from the background, the task provides options to extend or grow
+the object boundaries.  This is done by making multiple passes where
+pixels which have not been identified as object pixels but which neighbor
+object pixels are assigned to the object which they neighbor in any of
+the eight directions.  Each pass can be thought of as adding a ring
+of new pixels following the boundary of the object from the previous
+pass.
+
+When a non-object pixel neighbors two or more object pixels it is
+assigned to the object with the greater "flux".  The flux is the sum
+of the pixel value deviations from the background.
+
+The parameter \fIngrow\fR selects the maximum number of growing iterations.
+The parameter \fIagrow\fR selects the maximum fractional increase in
+the number  of original detected object pixels.  The number of pixels
+is called the "area"  of the object.  The growing of an object stops
+when either maximum is exceedd at the end of a growing iteration.
+.ih
+EXAMPLES
+1.  The following is a test example with default parameters that can be run
+by anyone.  An artificial galaxy field image is generated with the task
+\fBmkexample\fR (the \fBartdata\fR package is assumed to already be loaded)
+and a mask is created with \fBobjmasks\fR.  The image is displayed with
+the object mask overlayed in colors.
+
+.nf
+    np> mkexample galfield galfield
+    Creating example galfield in image galfield ...
+    np> objmasks omtype=color
+    List of images or MEF files: galfield
+    List of output object masks: gfmask
+    ACE:
+      Image: galfield - Example artificial galaxy field
+      Set sky and sigma:
+	Determine sky and sigma by surface fits:
+	  start line = 1, end line = 512, step = 51.1
+	  xorder = 2, yorder = 2, xterms = half
+	  hclip = 2., lclip = 3.
+	Determine sky and sigma by block statistics:
+	  Number of blocks: 5 5
+	  Number of pixels per block: 100 100
+	  Number of subblocks: 10 10
+	  Number of pixels per subblock: 50 50
+      Detect objects:
+	Convolution:
+	       1.      1.      1.
+	       1.      1.      1.
+	       1.      1.      1.
+	422 objects detected
+      Grow objects: ngrow = 2, agrow = 2.
+      Write object mask: gfmask[pl,type=mask]
+    np> display galfield 1
+    z1=371.5644 z2=455.8792
+    np> display galfield 2 overlay=gfmask[pl] ocolors="+203" 
+    z1=371.5644 z2=455.8792
+.fi
+
+2.  In the first example there was no input mask.  The next example
+creates a new object mask using the first object mask as an input
+"bad pixel mask".  While this is not the usual usage of the bad pixel
+mask it does illustrate an interesting option.  Note that the mask
+values in the input mask are mapped to an output value of 1 in the
+"colors" output.  In this example the output is forced to be a pl
+file by using the explicit extension.
+
+.nf
+    np> objmasks omtype=colors mask=gfmask[pl]
+    List of images or MEF files (galfield): 
+    List of output object masks (gfmask): gfmask1.pl
+    ACE:
+      Image: galfield - Example artificial galaxy field
+      Bad pixel mask: gfmask.pl
+      Set sky and sigma:
+	Determine sky and sigma by surface fits:
+	  start line = 1, end line = 512, step = 51.1
+	  xorder = 2, yorder = 2, xterms = half
+	  hclip = 2., lclip = 3.
+	Determine sky and sigma by block statistics:
+	  Number of blocks: 5 5
+	  Number of pixels per block: 100 100
+	  Number of subblocks: 10 10
+	  Number of pixels per subblock: 50 50
+      Detect objects:
+	Convolution:
+	       1.      1.      1.
+	       1.      1.      1.
+	       1.      1.      1.
+	44 objects detected
+      Grow objects: ngrow = 2, agrow = 2.
+      Write object mask: gfmask1.pl
+    np> display galfield 2 overlay=gfmask1 ocolors="+203" 
+    z1=371.5644 z2=455.8792
+.fi
+
+3.  The next example illustrates use with a multiextension file.  The
+example is two realizations of the galfield artificial data.
+
+.nf
+    np> mkexamples galfield mef.fits[im1]
+    Creating example galfield in image mef[im1] ...
+    np> mkexamples galfield mef[im2,append] oseed=2
+    Creating example galfield in image mef[im2,append] ...
+    np> objmasks
+    List of images or MEF files (galfield): mef
+    List of output object masks (gfmask1.pl): mefmask
+    ACE:
+      Image: mef[im1] - Example artificial galaxy field
+      Set sky and sigma:
+	Determine sky and sigma by surface fits:
+	  start line = 1, end line = 512, step = 51.1
+	  xorder = 2, yorder = 2, xterms = half
+	  hclip = 2., lclip = 3.
+	Determine sky and sigma by block statistics:
+	  Number of blocks: 5 5
+	  Number of pixels per block: 100 100
+	  Number of subblocks: 10 10
+	  Number of pixels per subblock: 50 50
+      Detect objects:
+	Convolution:
+	       1.      1.      1.
+	       1.      1.      1.
+	       1.      1.      1.
+	422 objects detected
+      Grow objects: ngrow = 2, agrow = 2.
+      Write object mask: mefmask[im1,append,type=mask]
+    ACE:
+      Image: mef[im2] - Example artificial galaxy field
+      Set sky and sigma:
+	Determine sky and sigma by surface fits:
+	  start line = 1, end line = 512, step = 51.1
+	  xorder = 2, yorder = 2, xterms = half
+	  hclip = 2., lclip = 3.
+	Determine sky and sigma by block statistics:
+	  Number of blocks: 5 5
+	  Number of pixels per block: 100 100
+	  Number of subblocks: 10 10
+	  Number of pixels per subblock: 50 50
+      Detect objects:
+	Convolution:
+	       1.      1.      1.
+	       1.      1.      1.
+	       1.      1.      1.
+	410 objects detected
+      Grow objects: ngrow = 2, agrow = 2.
+      Write object mask: mefmask[im2,append,type=mask]
+    np> display mef[im1] 1 over=mefmask[im1]
+    z1=371.5644 z2=455.8792
+    np> display mef[im2] 2 over=mefmask[im2]
+    z1=371.5666 z2=455.7844
+.fi
+
+4.  This example shows outputing the sky information.
+
+.nf
+    np> objmasks galfield gfmask2 sky=gfsky2
+    ACE:
+      Image: galfield - Example artificial galaxy field
+      Set sky and sigma:
+	Determine sky and sigma by surface fits:
+	  start line = 1, end line = 512, step = 51.1
+	  xorder = 2, yorder = 2, xterms = half
+	  hclip = 2., lclip = 3.
+	Determine sky and sigma by block statistics:
+	  Number of blocks: 5 5
+	  Number of pixels per block: 100 100
+	  Number of subblocks: 10 10
+	  Number of pixels per subblock: 50 50
+	Write sky map: gfsky2
+      Detect objects:
+	Convolution:
+	       1.      1.      1.
+	       1.      1.      1.
+	       1.      1.      1.
+	422 objects detected
+	Update sky map: gfsky2
+      Grow objects: ngrow = 2, agrow = 2.
+      Write object mask: gfmask2[pl,append,type=mask]
+    np> imstat gfsky2
+    #               IMAGE      NPIX      MEAN    STDDEV       MIN       MAX
+		   gfsky2        25     401.1    0.4397     400.3     401.9
+.fi
+
+5.  This examples shows specifying the sky information as constant values.
+In this case we already know that the artificial image has a
+constant background of 400 and a sigma of 10.
+
+.nf
+    np> objmasks galfield gfmask3 sky=400 sigma=10
+    ACE:
+      Image: galfield - Example artificial galaxy field
+      Set sky and sigma:
+	Use constant input sky: 400.
+	Use constant input sigma: 10.
+      Detect objects:
+	Convolution:
+	       1.      1.      1.
+	       1.      1.      1.
+	       1.      1.      1.
+	432 objects detected
+      Grow objects: ngrow = 2, agrow = 2.
+      Write object mask: gfmask3[pl,append,type=mask]
+.fi
+
+.ih
+REVISIONS
+.le
+.ih
+SEE ALSO
+.endhelp
+266c266
+< fit to produce a function for <S>.
+---
+> fit to produce a constant function for <S>.
+273c273,274
+< weights are used to fit a 2D plane for both <I> and <S>.
+---
+> weights are used to fit a 2D plane for both <I> and <S>.  The <S> surface
+> is a constant in order to avoid potential negative sigma values.
+