Initial commit

24 files changed, 3842 insertions, 0 deletions

diff --git a/noao/imred/crutil/src/Revisions b/noao/imred/crutil/src/Revisions
new file mode 100644
index 00000000..a5ef52f1
--- /dev/null
+++ b/noao/imred/crutil/src/Revisions
@@ -0,0 +1,151 @@
+.help revisions Nov99 crutil
+.nf
+t_cosmicrays.x
+    A pointer to a an array of pointers was used in one place as a real.  This
+    is an error when integer and real arrays are not of the same size; i.e.
+    on 64-bit architectures.  (8/2/12, Valdes)
+
+=======
+V2.16.1
+=======
+
+t_craverage.x
+    The pointer for the input mask buffer when using the same mask for both
+    input and output was not being read to include the edge lines which are
+    used during the computation resulting in a segmentation violation.
+    (11/24/04, Valdes)
+
+=======
+V2.12.2
+=======
+
+t_craverage.x
+    To grow an output mask create by the task requires closing the new mask
+    and reopening it READ_WRITE.  (10/23/02, Valdes)
+
+=======
+V2.12.1
+=======
+
+=====
+V2.12
+=====
+
+t_craverage.x
+t_crmedian.x
+    The mask name may be !<keyword>.  (3/22/02, Valdes)
+
+xtmaskname.x	+
+t_craverage.x
+t_crmedian.x
+t_crgrow.x
+mkpkg
+    Modified to allow FITS mask extensions.  (3/22/02, Valdes)
+    
+mkpkg
+    Added missing mkpkg depencies for several source files. (12/13/01, MJF)
+
+============================
+CRUTIL V1.5: August 22, 2001
+============================
+
+crcombine.cl
+crcombine.par
+../doc/crcombine.hlp
+    Modified to use new version of IMCOMBINE.  (8/22/01, Valdes)
+
+noao$imred/crutil/	+
+    Installed package into the NOAO.IMRED package.  (8/22/01, Valdes)
+
+t_craverage.x
+    The nrej features was done wrong for nrej of 1 or 2. (5/3/01, Valdes)
+
+t_craverage.x
+    1.  The amov and aclr calls were using the wrong buffer start.
+    2.  Added missing imtclose calls.
+    3.  The growing needed to be moved outside the block buffering.
+    (10/4/00, Valdes as diagnosed by Davis)
+
+t_craverage.x
+craverage.par
+../doc/craverage.hlp
+    Added an nrej parameter to excluded additional pixels from the average.
+    This is needed to deal with cosmic rays which are bigger than one
+    pixel or very nearby additional cosmic rays.  (9/13/00, Valdes)
+
+========================
+CRUTIL V1.4: Jan 6, 2000
+========================
+
+t_crmedian.x
+    The calculation of the sigma would reference uninitialized data if
+    the image size was not a multiple of the sigma block size.
+    (1/6/00, Valdes)
+
+t_craverage.x
+    The indexing of the pixels to use for the sigma calculation was wrong.
+    (1/6/00, Valdes)
+
+t_craverage.x	+
+craverage.par	+
+x_crutil.e
+mkpkg
+../crutil.cl
+../crutil.men
+../crutil.hd
+    New task that finds cosmic rays against an average excluding the candidate
+    pixel.  It also detects objects and prevents cosmic rays being detected
+    in them (i.e. the cores of stars).  (11/30/99, Valdes)
+
+t_crgrow.x
+crgrow.par
+    1.	Broke main grow loop into a subroutine that can be called by other
+	tasks.
+    2.  Added inval and outval parameters to allow selecting mask values.
+    (11/30/99, Valdes)
+
+========================
+CRUTIL V1.3: Oct 19, 1999
+========================
+
+crnebula.cl
+    The rin and rout parameters were not being used and instead
+    the default values were hardwired.  (10/19/99, Valdes)
+
+========================
+CRUTIL V1.2: Sep 4, 1998
+========================
+
+t_crmedian.x
+    On images with more than about 500K pixels the median operation is
+    done in overlapping blocks of lines.  The amount of overlap is
+    half of the line median size "lmed".  The bug is that on output
+    the overlap regions end up being zero.  The output i/o is now done
+    in non-overlapping blocks.  (9/4/98, Valdes)
+
+=====================
+CRUTIL V1.1: May 1998
+=====================
+
+crexamine.x
+cosmicrays.key	+
+cosmicrays.hlp
+    The graphical deletion and undeletion of candidates now includes the
+    keys 'e' and 'v' to delete and undelete from a marked rectangular
+    region.  Also the key file was moved to the source directory.
+    (Valdes/Shopbell 5/15/98)
+
+cosmicrays.par
+    Fixed type in "crmasks" parameter name.
+    (Valdes/Shopbell 5/15/98)
+
+===========
+CRUTIL V1.0
+===========
+    
+New package created April 24, 1998.
+
+=======
+V2.11.1
+=======
+.endhelp
diff --git a/noao/imred/crutil/src/cosmicrays.key b/noao/imred/crutil/src/cosmicrays.key
new file mode 100644
index 00000000..beac2835
--- /dev/null
+++ b/noao/imred/crutil/src/cosmicrays.key
@@ -0,0 +1,43 @@
+	   COSMIC RAY DETECTION AND REPLACEMENT
+
+   
+INTERACTIVE IMAGE CURSOR COMMANDS
+
+When using the image display cursor to define a list of training objects
+the following keystrokes may be given.
+
+?	Help
+c	Identify the object as a cosmic ray
+s	Identify the object as a star
+g	Switch to the graphics plot
+q	Quit and continue with the cleaning
+
+There are no colon commands.
+
+
+INTERACTIVE GRAPHICS CURSOR COMMANDS
+
+  The graph shows the ratio of the mean background subtracted flux within the
+detection window (excluding the candidate cosmic ray and the second brightest
+pixel) to the flux of the candidate cosmic ray pixel as a function of the
+mean flux.  Both coordinates have been multiplied by 100 so the ratio is in
+precent.  The peaks of extended objects have high flux ratios while true
+cosmic rays have low ratios.  The main purpose of this step is to set the
+flux ratio threshold for discriminating stars and galaxies.  The cursor
+keys are:
+
+?	Help
+a	Toggle between showing all objects and only the training objects
+d	Mark candidate for replacement (applys to '+' points)
+e	Mark candidates in a region for replacement (applys to '+' points)
+q	Quit.  Returns to training or to replace the selected pixels
+r	Redraw the graph
+s	Make a surface plot for the candidate nearest the cursor
+t	Set the flux ratio threshold at the y cursor position
+u	Mark candidate to not be replaced (applys to 'x' points)
+v	Mark candidates in a region to not be replaced (applys to 'x' points)
+w	Adjust the graph window (see GTOOLS)
+<space> Print the pixel coordinates for a candidate
+
+There are no colon commands except those for the windowing options (type
+:\help or see GTOOLS).
diff --git a/noao/imred/crutil/src/cosmicrays.par b/noao/imred/crutil/src/cosmicrays.par
new file mode 100644
index 00000000..9016d9f9
--- /dev/null
+++ b/noao/imred/crutil/src/cosmicrays.par
@@ -0,0 +1,17 @@
+input,s,a,,,,List of images in which to detect cosmic rays
+output,s,a,,,,List of cosmic ray replaced output images (optional)
+crmasks,s,h,"",,,"List of bad pixel masks (optional)
+"
+threshold,r,h,25.,,,Detection threshold above mean
+fluxratio,r,h,2.,,,Flux ratio threshold (in percent)
+npasses,i,h,5,1,,Number of detection passes
+window,s,h,"5","5|7",,"Size of detection window
+"
+interactive,b,h,yes,,,Examine parameters interactively?
+train,b,h,no,,,Use training objects?
+objects,*imcur,h,"",,,Cursor list of training objects
+savefile,f,h,"",,,File to save train objects
+plotfile,f,h,"",,,Plot file
+graphics,s,h,"stdgraph",,,Interactive graphics output device
+cursor,*gcur,h,"",,,Graphics cursor input
+answer,s,q,,"no|yes|NO|YES",,Review parameters for a particular image?
diff --git a/noao/imred/crutil/src/craverage.par b/noao/imred/crutil/src/craverage.par
new file mode 100644
index 00000000..b67c4585
--- /dev/null
+++ b/noao/imred/crutil/src/craverage.par
@@ -0,0 +1,23 @@
+input,f,a,,,,List of input images
+output,f,a,,,,List of output images
+crmask,f,h,,,,List of output cosmic ray and object masks
+average,f,h,,,,List of output block average filtered images
+sigma,f,h,,,,"List of output sigma images
+"
+navg,i,h,5,3,,Block average box size
+nrej,i,h,0,0,,Number of high pixels to reject from the average
+nbkg,i,h,5,1,,Background annulus width
+nsig,i,h,25,10,,Box size for sigma calculation
+var0,r,h,0.,0.,,Variance coefficient for DN^0 term
+var1,r,h,0.,0.,,Variance coefficient for DN^1 term
+var2,r,h,0.,0.,,"Variance coefficient for DN^2 term
+"
+crval,i,h,1,0,,Mask value for cosmic rays
+lcrsig,r,h,10.,,,Low cosmic ray sigma outside object
+hcrsig,r,h,5.,,,High cosmic ray sigma outside object
+crgrow,r,h,0.,0.,,"Cosmic ray grow radius
+"
+objval,i,h,0,0,,Mask value for objects
+lobjsig,r,h,10.,,,Low object detection sigma
+hobjsig,r,h,5.,,,High object detection sigma
+objgrow,r,h,0.,0.,,Object grow radius
diff --git a/noao/imred/crutil/src/crcombine.cl b/noao/imred/crutil/src/crcombine.cl
new file mode 100644
index 00000000..7c6bdc77
--- /dev/null
+++ b/noao/imred/crutil/src/crcombine.cl
@@ -0,0 +1,17 @@
+# CRCOMBINE -- Reject cosmic rays by combining multiple exposures.
+
+procedure crcombine (input, output)
+
+begin
+	imcombine (input, output, logfile=logfile, combine=combine,
+	    reject=reject, scale=scale, zero=zero, statsec=statsec,
+	    lsigma=lsigma, hsigma=hsigma, rdnoise=rdnoise, gain=gain,
+	    grow=grow, headers=headers, bpmasks=bpmasks, rejmasks=rejmasks,
+	    nrejmasks=nrejmasks, expmasks=expmasks, sigmas=sigmas,
+	    project=project, outtype=outtype, outlimits=outlimits,
+	    offsets=offsets, masktype=masktype, maskvalue=maskvalue,
+	    blank=blank, weight=weight, expname=expname,
+	    lthreshold=lthreshold, hthreshold=hthreshold, nlow=nlow,
+	    nhigh=nhigh, nkeep=nkeep, mclip=mclip, snoise=snoise,
+	    sigscale=sigscale, pclip=pclip)
+end
diff --git a/noao/imred/crutil/src/crcombine.par b/noao/imred/crutil/src/crcombine.par
new file mode 100644
index 00000000..95ae95ae
--- /dev/null
+++ b/noao/imred/crutil/src/crcombine.par
@@ -0,0 +1,45 @@
+# CRCOMBINE -- Image combine parameters
+
+input,s,a,,,,List of images to combine
+output,s,a,,,,List of output images
+logfile,s,h,"STDOUT",,,"Log file
+
+# Cosmic ray rejection parameters"
+combine,s,h,"average","average|median|sum",,Type of combine operation
+reject,s,h,"crreject","none|minmax|ccdclip|crreject|sigclip|avsigclip|pclip",,Type of rejection
+scale,s,h,"mode",,,Image scaling
+zero,s,h,"none",,,Image zero point offset
+statsec,s,h,"",,,Image section for computing statistics
+lsigma,r,h,10.,0.,,Lower sigma clipping factor
+hsigma,r,h,3.,0.,,Upper sigma clipping factor
+rdnoise,s,h,"0.",,,CCD readout noise (electrons)
+gain,s,h,"1.",,,CCD gain (electrons/DN)
+grow,r,h,0.,0.,,"Radius (pixels) for neighbor rejection
+
+# Additional output"
+headers,s,h,"",,,List of header files (optional)
+bpmasks,s,h,"",,,List of bad pixel masks (optional)
+rejmasks,s,h,"",,,List of rejection masks (optional)
+nrejmasks,s,h,"",,,List of number rejected masks (optional)
+expmasks,s,h,"",,,List of exposure masks (optional)
+sigmas,s,h,"",,,"List of sigma images (optional)
+
+# Additional parameters"
+project,b,h,no,,,Project highest dimension of input images?
+outtype,s,h,"real","short|ushort|integer|long|real|double",,Output image pixel datatype
+outlimits,s,h,"",,,Output limits (x1 x2 y1 y2 ...)
+offsets,f,h,"none",,,Input image offsets
+masktype,s,h,"none","none|goodvalue|badvalue|goodbits|badbits",,Mask type
+maskvalue,r,h,0,,,Mask value
+blank,r,h,0.,,,Value if there are no pixels
+weight,s,h,"none",,,Image weights
+expname,s,h,"",,,Image header exposure time keyword
+lthreshold,r,h,INDEF,,,Lower threshold
+hthreshold,r,h,INDEF,,,Upper threshold
+nlow,i,h,1,0,,minmax: Number of low pixels to reject
+nhigh,i,h,1,0,,minmax: Number of high pixels to reject
+nkeep,i,h,1,,,Minimum to keep (pos) or maximum to reject (neg)
+mclip,b,h,yes,,,Use median in sigma clipping algorithms?
+snoise,s,h,"0.",,,Sensitivity noise (fraction)
+sigscale,r,h,0.1,0.,,Tolerance for sigma clipping scaling corrections
+pclip,r,h,-0.5,,,pclip: Percentile clipping parameter
diff --git a/noao/imred/crutil/src/credit.cl b/noao/imred/crutil/src/credit.cl
new file mode 100644
index 00000000..8b844f33
--- /dev/null
+++ b/noao/imred/crutil/src/credit.cl
@@ -0,0 +1,13 @@
+# CREDIT -- Edit cosmic rays with an image display.
+
+procedure credit (input, output)
+
+begin
+	imedit (input, output, cursor=cursor, logfile=logfile,
+	    display=display, autodisplay=autodisplay,
+	    autosurface=autosurface, aperture=aperture, radius=radius,
+	    search=search, buffer=buffer, width=width, xorder=xorder,
+	    yorder=yorder, value=value, sigma=sigma, angh=angh, angv=angv,
+	    command=command, graphics=graphics, default=default,
+	    fixpix=fixpix)
+end
diff --git a/noao/imred/crutil/src/credit.par b/noao/imred/crutil/src/credit.par
new file mode 100644
index 00000000..a23e0d35
--- /dev/null
+++ b/noao/imred/crutil/src/credit.par
@@ -0,0 +1,22 @@
+input,s,a,,,,Images to be edited
+output,s,a,,,,Output images
+cursor,*imcur,h,"",,,Cursor input
+logfile,s,h,"",,,Logfile for record of cursor commands
+display,b,h,yes,,,Display images?
+autodisplay,b,h,yes,,,Automatic image display?
+autosurface,b,h,no,,,Automatic surface plots?
+aperture,s,h,"circular","|circular|square|",,Aperture type
+radius,r,h,2.,,,Substitution radius
+search,r,h,2.,,,Search radius
+buffer,r,h,1.,0.,,Background buffer width
+width,r,h,2.,1.,,Background width
+xorder,i,h,2,0,,Background x order
+yorder,i,h,2,0,,Background y order
+value,r,h,0.,,,Constant value substitution
+sigma,r,h,INDEF,,,Added noise sigma
+angh,r,h, -33.,,,Horizontal viewing angle (degrees)
+angv,r,h,25.,,,Vertical viewing angle (degrees)
+command,s,h,"display $image 1 erase=$erase fill=yes order=0 >& dev$null",,,Display command
+graphics,s,h,"stdgraph",,,Graphics device
+default,s,h,"b",,,Default option for x-y input
+fixpix,b,h,no,,,Fixpix style input?
diff --git a/noao/imred/crutil/src/crexamine.x b/noao/imred/crutil/src/crexamine.x
new file mode 100644
index 00000000..14e8d589
--- /dev/null
+++ b/noao/imred/crutil/src/crexamine.x
@@ -0,0 +1,626 @@
+include	<error.h>
+include	<syserr.h>
+include <imhdr.h>
+include	<gset.h>
+include	<mach.h>
+include	<pkg/gtools.h>
+include	"crlist.h"
+
+# CR_EXAMINE  -- Examine cosmic ray candidates interactively.
+# CR_GRAPH    -- Make a graph
+# CR_NEAREST  -- Find the nearest cosmic ray to the cursor.
+# CR_DELETE   -- Set replace flag for cosmic ray candidate nearest cursor.
+# CR_UNDELETE -- Set no replace flag for cosmic ray candidate nearest cursor.
+# CR_UPDATE   -- Change replacement flags, thresholds, and graphs.
+# CR_PLOT     -- Make log plot
+
+define	HELP	"crutil$src/cosmicrays.key"
+define	PROMPT	"cosmic ray options"
+
+# CR_EXAMINE -- Examine cosmic ray candidates interactively.
+
+procedure cr_examine (cr, gp, gt, im, fluxratio, first)
+
+pointer	cr			# Cosmic ray list
+pointer	gp			# GIO pointer
+pointer	gt			# GTOOLS pointer
+pointer	im			# Image pointer
+real	fluxratio		# Flux ratio threshold
+int	first			# Initial key
+
+char	cmd[SZ_LINE]
+int	i, newgraph, wcs, key, nc, nl, c1, c2, l1, l2, show
+real	wx, wy, x1, y1, x2, y2
+pointer	data
+
+int	clgcur()
+pointer	imgs2r()
+
+begin
+	# Set up the graphics.
+	call gt_sets (gt, GTPARAMS, IM_TITLE(im))
+
+	# Set image limits
+	nc = IM_LEN(im, 1)
+	nl = IM_LEN(im, 2)
+
+	# Enter cursor loop.
+	key = first
+	repeat {
+	    switch (key) {
+	    case '?': # Print help text.
+		call gpagefile (gp, HELP, PROMPT)
+	    case ':': # Colon commands.
+		switch (cmd[1]) {
+		case '/':
+		    call gt_colon (cmd, gp, gt, newgraph)
+		default:
+		    call printf ("\007")
+		}
+	    case 'a': # Toggle show all
+		if (show == 0)
+		    show = 1
+		else
+		    show = 0
+		newgraph = YES
+	    case 'd': # Delete candidate
+		call cr_delete (gp, wx, wy, cr, i, show)
+	    case 'e': # Delete candidates in region
+                x1 = wx; y1 = wy
+                call printf ("again:")
+                if (clgcur ("cursor", x2, y2, wcs, key, cmd, SZ_LINE) == EOF)
+                    return
+		call cr_delete_reg (gp, x1, y1, x2, y2, cr, show)
+	    case 'q': # Quit
+		break
+	    case 'r': # Redraw the graph.
+		newgraph = YES
+	    case 's': # Make surface plots
+		call cr_nearest (gp, wx, wy, cr, i, show)
+		c1 = max (1, int (Memr[CR_COL(cr)+i-1]) - 5)
+		c2 = min (nc, int (Memr[CR_COL(cr)+i-1]) + 5)
+		l1 = max (1, int (Memr[CR_LINE(cr)+i-1]) - 5)
+		l2 = min (nl, int (Memr[CR_LINE(cr)+i-1]) + 5)
+		data = imgs2r (im, c1, c2, l1, l2)
+		call gclear (gp)
+		call gsview (gp, 0.03, 0.48, 0.53, 0.98)
+		call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, -33., 25.)
+		call gsview (gp, 0.53, 0.98, 0.53, 0.98)
+		call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, -123., 25.)
+		call gsview (gp, 0.03, 0.48, 0.03, 0.48)
+		call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, 57., 25.)
+		call gsview (gp, 0.53, 0.98, 0.03, 0.48)
+		call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, 147., 25.)
+		call fprintf (STDERR, "[Type any key to continue]")
+		i = clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE)
+		newgraph = YES
+	    case 't': # Set threshold
+		call cr_update (gp, wy, cr, fluxratio, show)
+		call clputr ("fluxratio", fluxratio)
+	    case 'u': # Undelete candidate
+		call cr_undelete (gp, wx, wy, cr, i, show)
+	    case 'v': # Undelete candidates in region
+                x1 = wx; y1 = wy
+                call printf ("again:")
+                if (clgcur ("cursor", x2, y2, wcs, key, cmd, SZ_LINE) == EOF)
+                    return
+		call cr_undelete_reg (gp, x1, y1, x2, y2, cr, show)
+	    case 'w':# Window the graph.
+		call gt_window (gt, gp, "cursor", newgraph)
+	    case ' ': # Print info
+		call cr_nearest (gp, wx, wy, cr, i, show)
+		call printf ("%d %d\n")
+		    call pargr (Memr[CR_COL(cr)+i-1])
+		    call pargr (Memr[CR_LINE(cr)+i-1])
+	    case 'z': # NOP
+		newgraph = NO
+	    default: # Ring bell for unrecognized commands.
+		call printf ("\007")
+	    }
+
+	    # Update the graph if needed.
+	    if (newgraph == YES) {
+		call cr_graph (gp, gt, cr, fluxratio, show)
+	        newgraph = NO
+	    }
+	} until (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+end
+
+
+# CR_GRAPH -- Make a graph
+
+procedure cr_graph (gp, gt, cr, fluxratio, show)
+
+pointer	gp		# GIO pointer
+pointer	gt		# GTOOLS pointers
+pointer	cr		# Cosmic ray list
+real	fluxratio	# Flux ratio threshold
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr
+real	x1, x2, y1, y2
+pointer	sp, x, y, w, flag, index
+
+begin
+	call smark (sp)
+
+	call cr_show (show, cr, x, y, w, flag, index, ncr)
+	if (ncr == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	call gclear (gp)
+	call gt_ascale (gp, gt, Memr[x+1], Memr[y+1], ncr)
+	call gt_swind (gp, gt)
+	call gt_labax (gp, gt)
+
+	do i = 1, ncr {
+	    if ((Memi[flag+i] == NO) || (Memi[flag+i] == ALWAYSNO))
+		call gmark (gp, Memr[x+i], Memr[y+i], GM_PLUS, 2., 2.)
+	    else
+		call gmark (gp, Memr[x+i], Memr[y+i], GM_CROSS, 2., 2.)
+	    if (Memr[w+i] != 0.)
+		call gmark (gp, Memr[x+i], Memr[y+i], GM_BOX, 2., 2.)
+	}
+
+	call ggwind (gp, x1, x2, y1, y2)
+	call gseti (gp, G_PLTYPE, 2)
+	call gline (gp, x1, fluxratio, x2, fluxratio)
+
+	call sfree (sp)
+end
+
+
+# CR_NEAREST -- Find the nearest cosmic ray to the cursor.
+
+procedure cr_nearest (gp, wx, wy, cr, nearest, show)
+
+pointer	gp		# GIO pointer
+real	wx, wy		# Cursor position
+pointer	cr		# Cosmic ray list
+int	nearest		# Index of nearest point (returned)
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr
+real	x0, y0, x1, y1, x2, y2, r2, r2min
+pointer	sp, x, y, w, flag, index
+
+begin
+	call smark (sp)
+
+	call cr_show (show, cr, x, y, w, flag, index, ncr)
+	if (ncr == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Search for nearest point in NDC.
+	r2min = MAX_REAL
+	call gctran (gp, wx, wy, wx, wy, 1, 0)
+	do i = 1, ncr {
+	    x1 = Memr[x+i]
+	    y1 = Memr[y+i]
+	    call gctran (gp, x1, y1, x0, y0, 1, 0)
+	    r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	    if (r2 < r2min) {
+		r2min = r2
+		x2 = x1
+		y2 = y1
+		nearest = i
+	    }
+	}
+	if (index != NULL)
+	    nearest = Memi[index+nearest]
+
+	# Move the cursor to the selected point.
+	call gscur (gp, x2, y2)
+
+	call sfree (sp)
+end
+
+
+# CR_DELETE -- Set replace flag for cosmic ray candidate nearest cursor.
+
+procedure cr_delete (gp, wx, wy, cr, nearest, show)
+
+pointer	gp		# GIO pointer
+real	wx, wy		# Cursor position
+pointer	cr		# Cosmic ray list
+int	nearest		# Index of nearest point (returned)
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr
+real	x0, y0, x1, y1, x2, y2, r2, r2min
+pointer	sp, x, y, w, flag, index
+
+begin
+	call smark (sp)
+
+	call cr_show (show, cr, x, y, w, flag, index, ncr)
+	if (ncr == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Search for nearest point in NDC.
+	nearest = 0
+	r2min = MAX_REAL
+	call gctran (gp, wx, wy, wx, wy, 1, 0)
+	do i = 1, ncr {
+	    if ((Memi[flag+i] == YES) || (Memi[flag+i] == ALWAYSYES))
+		next
+	    x1 = Memr[x+i]
+	    y1 = Memr[y+i]
+	    call gctran (gp, x1, y1, x0, y0, 1, 0)
+	    r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	    if (r2 < r2min) {
+		r2min = r2
+		x2 = x1
+		y2 = y1
+		nearest = i
+	    }
+	}
+
+	# Move the cursor to the selected point and mark the deleted point.
+	if (nearest > 0) {
+	    if (index != NULL)
+		nearest = Memi[index+nearest]
+	    Memi[CR_FLAG(cr)+nearest-1] = ALWAYSYES
+	    Memi[CR_WT(cr)+nearest-1] = -1
+	    call gscur (gp, x2, y2)
+	    call gseti (gp, G_PMLTYPE, 0)
+	    y2 = Memr[CR_RATIO(cr)+nearest-1]
+	    call gmark (gp, x2, y2, GM_PLUS, 2., 2.)
+	    call gseti (gp, G_PMLTYPE, 1)
+	    call gmark (gp, x2, y2, GM_CROSS, 2., 2.)
+	}
+
+	call sfree (sp)
+end
+
+
+# CR_DELETE_REG -- Set replace flag for cosmic ray candidates in a region.
+
+procedure cr_delete_reg (gp, wx1, wy1, wx2, wy2, cr, show)
+
+pointer	gp			# GIO pointer
+real	wx1, wy1, wx2, wy2	# Cursor positions
+pointer	cr			# Cosmic ray list
+int	show			# Show (0=all, 1=train)
+
+int	i, j, ncr
+real	x0, y0, x1, y1
+pointer	sp, x, y, w, flag, index
+
+begin
+	call smark (sp)
+
+	call cr_show (show, cr, x, y, w, flag, index, ncr)
+	if (ncr == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Check order of region points.
+	if (wx1 > wx2) {
+	    x0 = wx1
+	    wx1 = wx2
+	    wx2 = x0
+	}
+	if (wy1 > wy2) {
+	    y0 = wy1
+	    wy1 = wy2
+	    wy2 = y0
+	}
+
+	# Check if point in region.
+	call gctran (gp, wx1, wy1, wx1, wy1, 1, 0)
+	call gctran (gp, wx2, wy2, wx2, wy2, 1, 0)
+	do i = 1, ncr {
+	    if ((Memi[flag+i] == YES) || (Memi[flag+i] == ALWAYSYES))
+		next
+	    x1 = Memr[x+i]
+	    y1 = Memr[y+i]
+	    call gctran (gp, x1, y1, x0, y0, 1, 0)
+
+	    # Mark the deleted points.
+	    if ((x0 > wx1) && (x0 < wx2) && (y0 > wy1) && (y0 < wy2)) {
+		if (index != NULL)
+		    j = Memi[index+i]
+		else
+		    j = i
+		Memi[CR_FLAG(cr)+j-1] = ALWAYSYES
+		Memi[CR_WT(cr)+j-1] = -1
+		call gscur (gp, x1, y1)
+		call gseti (gp, G_PMLTYPE, 0)
+		y1 = Memr[CR_RATIO(cr)+j-1]
+		call gmark (gp, x1, y1, GM_PLUS, 2., 2.)
+		call gseti (gp, G_PMLTYPE, 1)
+		call gmark (gp, x1, y1, GM_CROSS, 2., 2.)
+	    }
+	}
+	call sfree (sp)
+end
+
+
+# CR_UNDELETE -- Set no replace flag for cosmic ray candidate nearest cursor.
+
+procedure cr_undelete (gp, wx, wy, cr, nearest, show)
+
+pointer	gp		# GIO pointer
+real	wx, wy		# Cursor position
+pointer	cr		# Cosmic ray list
+int	nearest		# Index of nearest point (returned)
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr
+real	x0, y0, x1, y1, x2, y2, r2, r2min
+pointer	sp, x, y, w, flag, index
+
+begin
+	call smark (sp)
+
+	call cr_show (show, cr, x, y, w, flag, index, ncr)
+	if (ncr == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Search for nearest point in NDC.
+	nearest = 0
+	r2min = MAX_REAL
+	call gctran (gp, wx, wy, wx, wy, 1, 0)
+	do i = 1, ncr {
+	    if ((Memi[flag+i] == NO) || (Memi[flag+i] == ALWAYSNO))
+		next
+	    x1 = Memr[x+i]
+	    y1 = Memr[y+i]
+	    call gctran (gp, x1, y1, x0, y0, 1, 0)
+	    r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	    if (r2 < r2min) {
+		r2min = r2
+		x2 = x1
+		y2 = y1
+		nearest = i
+	    }
+	}
+
+	# Move the cursor to the selected point and mark the delete point.
+	if (nearest > 0) {
+	    if (index != NULL)
+		nearest = Memi[index+nearest]
+	    Memi[CR_FLAG(cr)+nearest-1] = ALWAYSNO
+	    Memi[CR_WT(cr)+nearest-1] = 1
+	    call gscur (gp, x2, y2)
+
+	    call gseti (gp, G_PMLTYPE, 0)
+	    y2 = Memr[CR_RATIO(cr)+nearest-1]
+	    call gmark (gp, x2, y2, GM_CROSS, 2., 2.)
+	    call gseti (gp, G_PMLTYPE, 1)
+	    call gmark (gp, x2, y2, GM_PLUS, 2., 2.)
+	}
+
+	call sfree (sp)
+end
+
+
+# CR_UNDELETE_REG -- Set no replace flag for cosmic ray candidates in a region.
+
+procedure cr_undelete_reg (gp, wx1, wy1, wx2, wy2, cr, show)
+
+pointer	gp			# GIO pointer
+real	wx1, wy1, wx2, wy2	# Cursor positions
+pointer	cr			# Cosmic ray list
+int	show			# Show (0=all, 1=train)
+
+int	i, j, ncr
+real	x0, y0, x1, y1
+pointer	sp, x, y, w, flag, index
+
+begin
+	call smark (sp)
+
+	call cr_show (show, cr, x, y, w, flag, index, ncr)
+	if (ncr == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Check order of region points.
+	if (wx1 > wx2) {
+	    x0 = wx1
+	    wx1 = wx2
+	    wx2 = x0
+	}
+	if (wy1 > wy2) {
+	    y0 = wy1
+	    wy1 = wy2
+	    wy2 = y0
+	}
+
+	# Check if point in region.
+	call gctran (gp, wx1, wy1, wx1, wy1, 1, 0)
+	call gctran (gp, wx2, wy2, wx2, wy2, 1, 0)
+	do i = 1, ncr {
+	    if ((Memi[flag+i] == NO) || (Memi[flag+i] == ALWAYSNO))
+		next
+	    x1 = Memr[x+i]
+	    y1 = Memr[y+i]
+	    call gctran (gp, x1, y1, x0, y0, 1, 0)
+
+	    # Mark the deleted points.
+	    if ((x0 > wx1) && (x0 < wx2) && (y0 > wy1) && (y0 < wy2)) {
+		if (index != NULL)
+		    j = Memi[index+i]
+		else
+		    j = i
+		Memi[CR_FLAG(cr)+j-1] = ALWAYSNO
+		Memi[CR_WT(cr)+j-1] = 1
+		call gscur (gp, x1, y1)
+		call gseti (gp, G_PMLTYPE, 0)
+		y1 = Memr[CR_RATIO(cr)+j-1]
+		call gmark (gp, x1, y1, GM_CROSS, 2., 2.)
+		call gseti (gp, G_PMLTYPE, 1)
+		call gmark (gp, x1, y1, GM_PLUS, 2., 2.)
+	    }
+	}
+	call sfree (sp)
+end
+
+
+# CR_UPDATE -- Change replacement flags, thresholds, and graphs.
+
+procedure cr_update (gp, wy, cr, fluxratio, show)
+
+pointer	gp		# GIO pointer
+real	wy		# Y cursor position
+pointer	cr		# Cosmic ray list
+real	fluxratio	# Flux ratio threshold
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr, flag
+real	x1, x2, y1, y2
+pointer	x, y, f
+
+begin
+	call gseti (gp, G_PLTYPE, 0)
+	call ggwind (gp, x1, x2, y1, y2)
+	call gline (gp, x1, fluxratio, x2, fluxratio)
+	fluxratio = wy
+	call gseti (gp, G_PLTYPE, 2)
+	call gline (gp, x1, fluxratio, x2, fluxratio)
+
+	if (show == 1)
+	    return
+
+	ncr = CR_NCR(cr)
+	x = CR_FLUX(cr) - 1
+	y = CR_RATIO(cr) - 1
+	f = CR_FLAG(cr) - 1
+
+	do i = 1, ncr {
+	    flag = Memi[f+i]
+	    if ((flag == ALWAYSYES) || (flag == ALWAYSNO))
+		next
+	    x1 = Memr[x+i]
+	    y1 = Memr[y+i]
+	    if (flag == NO) {
+	        if (y1 < fluxratio) {
+		    Memi[f+i] = YES
+		    call gseti (gp, G_PMLTYPE, 0)
+		    call gmark (gp, x1, y1, GM_PLUS, 2., 2.)
+		    call gseti (gp, G_PMLTYPE, 1)
+		    call gmark (gp, x1, y1, GM_CROSS, 2., 2.)
+		}
+	    } else {
+	        if (y1 >= fluxratio) {
+		    Memi[f+i] = NO
+		    call gseti (gp, G_PMLTYPE, 0)
+		    call gmark (gp, x1, y1, GM_CROSS, 2., 2.)
+		    call gseti (gp, G_PMLTYPE, 1)
+		    call gmark (gp, x1, y1, GM_PLUS, 2., 2.)
+		}
+	    }
+	}
+end
+
+
+# CR_PLOT -- Make log plot
+
+procedure cr_plot (cr, im, fluxratio)
+
+pointer	cr			# Cosmic ray list
+pointer	im			# Image pointer
+real	fluxratio		# Flux ratio threshold
+
+int	fd, open(), errcode()
+pointer	sp, fname, gp, gt, gopen(), gt_init()
+errchk	gopen
+
+begin
+	call smark (sp)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+
+	# Open the plotfile.
+	call clgstr ("plotfile", Memc[fname], SZ_FNAME)
+	iferr (fd = open (Memc[fname], APPEND, BINARY_FILE)) {
+	    if (errcode() != SYS_FNOFNAME)
+		call erract (EA_WARN)
+	    return
+	}
+
+	# Set up the graphics.
+	gp = gopen ("stdplot", NEW_FILE, fd)
+	gt = gt_init()
+	call gt_sets (gt, GTTYPE, "mark")
+	call gt_sets (gt, GTXTRAN, "log")
+	call gt_setr (gt, GTXMIN, 10.)
+	call gt_setr (gt, GTYMIN, 0.)
+	call gt_sets (gt, GTTITLE, "Parameters of cosmic rays candidates")
+	call gt_sets (gt, GTPARAMS, IM_TITLE(im))
+	call gt_sets (gt, GTXLABEL, "Flux")
+	call gt_sets (gt, GTYLABEL, "Flux Ratio")
+
+	call cr_graph (gp, gt, cr, fluxratio, 'r')
+
+	call gt_free (gt)
+	call gclose (gp)
+	call close (fd)
+	call sfree (sp)
+end
+
+
+# CR_SHOW -- Select data to show.
+# This returns pointers to the data.  Note the pointers are salloc from
+# the last smark which is done by the calling program.
+
+procedure cr_show (show, cr, x, y, w, flag, index, ncr)
+
+int	show		#I Data to show (0=all, 1=train)
+pointer	cr		#I CR data
+pointer	x		#O Fluxes
+pointer	y		#O Ratios
+pointer	w		#O Weights
+pointer	flag		#O Flags
+pointer	index		#O Index into CR data (if not null)
+int	ncr		#O Number of selected data points
+
+int	i
+
+begin
+	switch (show) {
+	case 0:
+	    ncr = CR_NCR(cr)
+	    x = CR_FLUX(cr) - 1
+	    y = CR_RATIO(cr) - 1
+	    w = CR_WT(cr) - 1
+	    flag = CR_FLAG(cr) - 1
+	    index = NULL
+	case 1:
+	    ncr = CR_NCR(cr)
+	    call salloc (x, ncr, TY_REAL)
+	    call salloc (y, ncr, TY_REAL)
+	    call salloc (w, ncr, TY_REAL)
+	    call salloc (flag, ncr, TY_INT)
+	    call salloc (index, ncr, TY_INT)
+
+	    ncr = 0
+	    x = x - 1
+	    y = y - 1
+	    w = w - 1
+	    flag = flag - 1
+	    index = index - 1
+
+	    do i = 1, CR_NCR(cr) {
+		if (Memr[CR_WT(cr)+i-1] == 0.)
+		    next
+		ncr = ncr + 1
+		Memr[x+ncr] = Memr[CR_FLUX(cr)+i-1]
+		Memr[y+ncr] = Memr[CR_RATIO(cr)+i-1]
+		Memr[w+ncr] = Memr[CR_WT(cr)+i-1]
+		Memi[flag+ncr] = Memi[CR_FLAG(cr)+i-1]
+		Memi[index+ncr] = i
+	    }
+	}
+end
diff --git a/noao/imred/crutil/src/crfind.x b/noao/imred/crutil/src/crfind.x
new file mode 100644
index 00000000..58850940
--- /dev/null
+++ b/noao/imred/crutil/src/crfind.x
@@ -0,0 +1,305 @@
+include	<math/gsurfit.h>
+
+# CR_FIND -- Find cosmic ray candidates.
+# This procedure is an interface to special procedures specific to a given
+# window size.
+
+procedure cr_find (cr, threshold, data, nc, nl, col, line,
+    sf1, sf2, x, y, z, w)
+
+pointer	cr						# Cosmic ray list
+real	threshold					# Detection threshold
+pointer	data[ARB]					# Data lines
+int	nc						# Number of columns
+int	nl						# Number of lines
+int	col						# First column
+int	line						# Center line
+pointer	sf1, sf2					# Surface fitting
+real	x[ARB], y[ARB], z[ARB], w[ARB]			# Surface arrays
+
+pointer	a, b, c, d, e, f, g
+
+begin
+	switch (nl) {
+	case 5:
+	    a = data[1]
+	    b = data[2]
+	    c = data[3]
+	    d = data[4]
+	    e = data[5]
+	    call cr_find5 (cr, threshold,  col, line, Memr[a], Memr[b],
+		Memr[c], Memr[d], Memr[e], nc, sf1, sf2, x, y, z, w)
+	case 7:
+	    a = data[1]
+	    b = data[2]
+	    c = data[3]
+	    d = data[4]
+	    e = data[5]
+	    f = data[6]
+	    g = data[7]
+	    call cr_find7 (cr, threshold,  col, line, Memr[a], Memr[b],
+		Memr[c], Memr[d], Memr[e], Memr[f], Memr[g], nc,
+		sf1, sf2, x, y, z, w)
+	}
+end
+
+
+# CR_FIND7 -- Find cosmic rays candidates in 7x7 window.
+# This routine finds cosmic rays candidates with the following algorithm.
+#    1. If the pixel is not a local maximum relative to it's 48 neighbors
+#	go on to the next pixel.
+#    2. Identify the next strongest pixel in the 7x7 region.
+#	This suspect pixel is excluded in the following.
+#    2. Compute the flux of the 7x7 region excluding the cosmic ray
+#	candidate and the suspect pixel.
+#    3. The candidate must exceed the average flux per pixel by a specified
+#	threshold.  If not go on to the next pixel.
+#    4. Fit a plane to the border pixels (excluding the suspect pixel).
+#    5. Subtract the background defined by the plane.
+#    6. Determine a replacement value as the average of the four adjacent
+#	pixels (excluding the suspect pixels).
+#    7. Add the pixel to the cosmic ray candidate list.
+
+procedure cr_find7 (cr, threshold, col, line, a, b, c, d, e, f, g, n,
+	sf1, sf2, x, y, z, w)
+
+pointer	cr						# Cosmic ray list
+real	threshold					# Detection threshold
+int	col						# First column
+int	line						# Line
+real	a[ARB], b[ARB], c[ARB], d[ARB]			# Image lines
+real	e[ARB], f[ARB], g[ARB]				# Image lines
+int	n						# Number of columns
+pointer	sf1, sf2					# Surface fitting
+real	x[49], y[49], z[49], w[49]			# Surface arrays
+
+real	bkgd[49] 
+int	i1, i2, i3, i4, i5, i6, i7, j, j1, j2
+real	p, flux, replace, asumr()
+pointer	sf
+
+begin
+	for (i4=4; i4<=n-3; i4=i4+1) {
+	    # Must be local maxima.
+	    p = d[i4]
+	    if (p<a[i4]||p<b[i4]||p<c[i4]||p<e[i4]||p<f[i4]||p<g[i4])
+		next
+	    i1 = i4 - 3
+	    if (p<a[i1]||p<b[i1]||p<c[i1]||p<d[i1]||p<e[i1]||p<f[i1]||p<g[i1])
+		next
+	    i2 = i4 - 2
+	    if (p<a[i2]||p<b[i2]||p<c[i2]||p<d[i2]||p<e[i2]||p<f[i2]||p<g[i2])
+		next
+	    i3 = i4 - 1
+	    if (p<a[i3]||p<b[i3]||p<c[i3]||p<d[i3]||p<e[i3]||p<f[i3]||p<g[i3])
+		next
+	    i5 = i4 + 1
+	    if (p<a[i5]||p<b[i5]||p<c[i5]||p<d[i5]||p<e[i5]||p<f[i5]||p<g[i5])
+		next
+	    i6 = i4 + 2
+	    if (p<a[i6]||p<b[i6]||p<c[i6]||p<d[i6]||p<e[i6]||p<f[i6]||p<g[i6])
+		next
+	    i7 = i4 + 3
+	    if (p<a[i7]||p<b[i7]||p<c[i7]||p<d[i7]||p<e[i7]||p<f[i7]||p<g[i7])
+		next
+
+	    # Convert to a single array in surface fitting order. 
+	    call amovr (a[i1], z[1], 7)
+	    z[8] = b[i7]; z[9] = c[i7]; z[10] = d[i7]; z[11] = e[i7]
+	    z[12] = f[i7]; z[13] = g[i7]; z[14] = g[i6]; z[15] = g[i5]
+	    z[16] = f[i4]; z[17] = g[i3]; z[18] = g[i2]; z[19] = g[i1]
+	    z[20] = f[i1]; z[21] = e[i1]; z[22] = d[i1]; z[23] = c[i1]
+	    z[24] = b[i1]
+	    call amovr (b[i2], z[25], 5)
+	    call amovr (c[i2], z[30], 5)
+	    call amovr (d[i2], z[35], 5)
+	    call amovr (e[i2], z[40], 5)
+	    call amovr (f[i2], z[45], 5)
+
+	    # Find the highest point excluding the center.
+	    j1 = 37; j2 = 1
+	    do j = 2, 49 {
+		if (j == j1)
+		    next
+		if (z[j] > z[j2])
+		    j2 = j 
+	    }
+
+	    # Compute the flux excluding the extreme points.
+	    flux = (asumr (z, 49) - z[j1] - z[j2]) / 47
+
+	    # Pixel must be exceed specified threshold.
+	    if (p < flux + threshold)
+		next
+
+	    # Fit and subtract the background.
+	    if (j2 < 25) {
+	        w[j2] = 0
+		sf = sf2
+	        call gsfit (sf, x, y, z, w, 24, WTS_USER, j)
+		w[j2] = 1
+	    } else {
+		sf = sf1
+	        call gsrefit (sf, x, y, z, w, j)
+	    }
+
+	    call gsvector (sf, x, y, bkgd, 49)
+	    call asubr (z, bkgd, z, 49)
+	    p = z[j1]
+
+	    # Compute the flux excluding the extreme points.
+	    flux = (asumr (z, 49) - z[j1] - z[j2]) / 47
+
+	    # Determine replacement value from four nearest neighbors again
+	    # excluding the most deviant pixels.
+	    replace = 0
+	    j = 0
+	    if (j2 != 32) {
+		replace = replace + c[i4]
+		j = j + 1
+	    }
+	    if (j2 != 36) {
+		replace = replace + d[i3]
+		j = j + 1
+	    }
+	    if (j2 != 38) {
+		replace = replace + d[i5]
+		j = j + 1
+	    }
+	    if (j2 != 42) {
+		replace = replace + e[i4]
+		j = j + 1
+	    }
+	    replace = replace / j
+
+	    # Add pixel to cosmic ray list.
+	    flux = 100. * flux
+	    call cr_add (cr, col+i4-1, line, flux, flux/p, 0., replace, 0)
+	    i4 = i7
+	}
+end
+
+
+# CR_FIND5 -- Find cosmic rays candidates in 5x5 window.
+# This routine finds cosmic rays candidates with the following algorithm.
+#    1. If the pixel is not a local maximum relative to it's 24 neighbors
+#	go on to the next pixel.
+#    2. Identify the next strongest pixel in the 5x5 region.
+#	This suspect pixel is excluded in the following.
+#    2. Compute the flux of the 5x5 region excluding the cosmic ray
+#	candidate and the suspect pixel.
+#    3. The candidate must exceed the average flux per pixel by a specified
+#	threshold.  If not go on to the next pixel.
+#    4. Fit a plane to the border pixels (excluding the suspect pixel).
+#    5. Subtract the background defined by the plane.
+#    6. Determine a replacement value as the average of the four adjacent
+#	pixels (excluding the suspect pixels).
+#    7. Add the pixel to the cosmic ray candidate list.
+
+procedure cr_find5 (cr, threshold, col, line, a, b, c, d, e, n,
+	sf1, sf2, x, y, z, w)
+
+pointer	cr						# Cosmic ray list
+real	threshold					# Detection threshold
+int	col						# First column
+int	line						# Line
+real	a[ARB], b[ARB], c[ARB], d[ARB], e[ARB]		# Image lines
+int	n						# Number of columns
+pointer	sf1, sf2					# Surface fitting
+real	x[25], y[25], z[25], w[25]			# Surface arrays
+
+real	bkgd[25] 
+int	i1, i2, i3, i4, i5, j, j1, j2
+real	p, flux, replace, asumr()
+pointer	sf
+
+begin
+	for (i3=3; i3<=n-2; i3=i3+1) {
+	    # Must be local maxima.
+	    p = c[i3]
+	    if (p<a[i3]||p<b[i3]||p<d[i3]||p<e[i3])
+		next
+	    i1 = i3 - 2
+	    if (p<a[i1]||p<b[i1]||p<c[i1]||p<d[i1]||p<e[i1])
+		next
+	    i2 = i3 - 1
+	    if (p<a[i2]||p<b[i2]||p<c[i2]||p<d[i2]||p<e[i2])
+		next
+	    i4 = i3 + 1
+	    if (p<a[i4]||p<b[i4]||p<c[i4]||p<d[i4]||p<e[i4])
+		next
+	    i5 = i3 + 2
+	    if (p<a[i5]||p<b[i5]||p<c[i5]||p<d[i5]||p<e[i5])
+		next
+
+	    # Convert to a single array in surface fitting order. 
+	    call amovr (a[i1], z[1], 5)
+	    z[6] = b[i5]; z[7] = c[i5]; z[8] = d[i5]; z[9] = e[i5]
+	    z[10] = e[i4]; z[11] = e[i3]; z[12] = e[i2]; z[13] = e[i1]
+	    z[14] = d[i1]; z[15] = c[i1]; z[16] = b[i1]
+	    call amovr (b[i2], z[17], 3)
+	    call amovr (c[i2], z[20], 3)
+	    call amovr (d[i2], z[23], 3)
+
+	    # Find the highest point excluding the center.
+	    j1 = 21; j2 = 1
+	    do j = 2, 25 {
+		if (j == j1)
+		    next
+		if (z[j] > z[j2])
+		    j2 = j 
+	    }
+
+	    # Compute the flux excluding the extreme points.
+	    flux = (asumr (z, 25) - z[j1] - z[j2]) / 23
+
+	    # Pixel must be exceed specified threshold.
+	    if (p < flux + threshold)
+		next
+
+	    # Fit and subtract the background.
+	    if (j2 < 17) {
+	        w[j2] = 0
+		sf = sf2
+	        call gsfit (sf, x, y, z, w, 16, WTS_USER, j)
+		w[j2] = 1
+	    } else {
+		sf = sf1
+	        call gsrefit (sf, x, y, z, w, j)
+	    }
+
+	    call gsvector (sf, x, y, bkgd, 25)
+	    call asubr (z, bkgd, z, 25)
+	    p = z[j1]
+
+	    # Compute the flux excluding the extreme points.
+	    flux = (asumr (z, 25) - z[j1] - z[j2]) / 23
+
+	    # Determine replacement value from four nearest neighbors again
+	    # excluding the most deviant pixels.
+	    replace = 0
+	    j = 0
+	    if (j2 != 18) {
+		replace = replace + b[i3]
+		j = j + 1
+	    }
+	    if (j2 != 20) {
+		replace = replace + c[i2]
+		j = j + 1
+	    }
+	    if (j2 != 22) {
+		replace = replace + c[i4]
+		j = j + 1
+	    }
+	    if (j2 != 24) {
+		replace = replace + d[i3]
+		j = j + 1
+	    }
+	    replace = replace / j
+
+	    # Add pixel to cosmic ray list.
+	    flux = 100. * flux
+	    call cr_add (cr, col+i3-1, line, flux, flux/p, 0., replace, 0)
+	    i3 = i5
+	}
+end
diff --git a/noao/imred/crutil/src/crfix.cl b/noao/imred/crutil/src/crfix.cl
new file mode 100644
index 00000000..68947c7a
--- /dev/null
+++ b/noao/imred/crutil/src/crfix.cl
@@ -0,0 +1,20 @@
+# CRFIX -- Replace cosmic rays in an image using a cosmic ray mask.
+
+procedure crfix (input, output, crmask)
+
+file	input 			{prompt="Input image"}
+file	output 			{prompt="Output image"}
+file	crmask 			{prompt="Cosmic ray mask"}
+
+begin
+	file	in, out, crm
+
+	in = input
+	out = output
+	crm = crmask
+
+	if (in != out)
+	    imcopy (in, out, verbose=no)
+	fixpix (out, crm, linterp="INDEF", cinterp="INDEF", verbose=no,
+	    pixels=no)
+end
diff --git a/noao/imred/crutil/src/crgrow.par b/noao/imred/crutil/src/crgrow.par
new file mode 100644
index 00000000..f57eb3cc
--- /dev/null
+++ b/noao/imred/crutil/src/crgrow.par
@@ -0,0 +1,7 @@
+# CRGROW
+
+input,s,a,,,,Input cosmic ray masks
+output,s,a,,,,Output cosmic ray masks
+radius,r,h,1.,,,Grow radius
+inval,i,h,INDEF,,,Input mask value to grow (INDEF for any)
+outval,i,h,INDEF,,,Output grown mask value (INDEF for any)
diff --git a/noao/imred/crutil/src/crlist.h b/noao/imred/crutil/src/crlist.h
new file mode 100644
index 00000000..1ed498a7
--- /dev/null
+++ b/noao/imred/crutil/src/crlist.h
@@ -0,0 +1,17 @@
+define	CR_ALLOC	100		# Allocation block size
+define	CR_LENSTRUCT	9		# Length of structure
+
+define	CR_NCR		Memi[$1]	# Number of cosmic rays
+define	CR_NALLOC	Memi[$1+1]	# Length of cosmic ray list
+define	CR_COL		Memi[$1+2]	# Pointer to columns
+define	CR_LINE		Memi[$1+3]	# Pointer to lines
+define	CR_FLUX		Memi[$1+4]	# Pointer to fluxes
+define	CR_RATIO	Memi[$1+5]	# Pointer to flux ratios
+define	CR_WT		Memi[$1+6]	# Pointer to training weights
+define	CR_REPLACE	Memi[$1+7]	# Pointer to replacement values
+define	CR_FLAG		Memi[$1+8]	# Pointer to rejection flag
+
+define	ALWAYSNO	3
+define	ALWAYSYES	4
+
+define	CR_RMAX		3.		# Maximum radius for matching
diff --git a/noao/imred/crutil/src/crlist.x b/noao/imred/crutil/src/crlist.x
new file mode 100644
index 00000000..bb49fb03
--- /dev/null
+++ b/noao/imred/crutil/src/crlist.x
@@ -0,0 +1,417 @@
+include	<error.h>
+include	<syserr.h>
+include	<gset.h>
+include	<pmset.h>
+include	"crlist.h"
+
+define	HELP	"noao$lib/scr/cosmicrays.key"
+define	PROMPT	"cosmic ray options"
+
+# CR_OPEN    -- Open cosmic ray list
+# CR_CLOSE   -- Close cosmic ray list
+# CR_ADD     -- Add a cosmic ray candidate to cosmic ray list.
+# CR_TRAIN   -- Set flux ratio threshold from a training set.
+# CR_FINDTHRESH -- Find flux ratio.
+# CR_WEIGHT  -- Compute the training weight at a particular flux ratio.
+# CR_FLAGS   -- Set cosmic ray reject flags.
+# CR_BADPIX  -- Store cosmic rays in bad pixel list.
+# CR_REPLACE -- Replace cosmic rays in image with replacement values.
+
+# CR_OPEN -- Open cosmic ray list
+
+procedure cr_open (cr)
+
+pointer	cr		# Cosmic ray list pointer
+errchk	malloc
+
+begin
+	call malloc (cr, CR_LENSTRUCT, TY_STRUCT)
+	call malloc (CR_COL(cr), CR_ALLOC, TY_REAL)
+	call malloc (CR_LINE(cr), CR_ALLOC, TY_REAL)
+	call malloc (CR_FLUX(cr), CR_ALLOC, TY_REAL)
+	call malloc (CR_RATIO(cr), CR_ALLOC, TY_REAL)
+	call malloc (CR_WT(cr), CR_ALLOC, TY_REAL)
+	call malloc (CR_REPLACE(cr), CR_ALLOC, TY_REAL)
+	call malloc (CR_FLAG(cr), CR_ALLOC, TY_INT)
+	CR_NCR(cr) = 0
+	CR_NALLOC(cr) = CR_ALLOC
+end
+
+
+# CR_CLOSE -- Close cosmic ray list
+
+procedure cr_close (cr)
+
+pointer	cr		# Cosmic ray list pointer
+
+begin
+	call mfree (CR_COL(cr), TY_REAL)
+	call mfree (CR_LINE(cr), TY_REAL)
+	call mfree (CR_FLUX(cr), TY_REAL)
+	call mfree (CR_RATIO(cr), TY_REAL)
+	call mfree (CR_WT(cr), TY_REAL)
+	call mfree (CR_REPLACE(cr), TY_REAL)
+	call mfree (CR_FLAG(cr), TY_INT)
+	call mfree (cr, TY_STRUCT)
+end
+
+# CR_ADD -- Add a cosmic ray candidate to cosmic ray list.
+
+procedure cr_add (cr, col, line, flux, ratio, wt, replace, flag)
+
+pointer	cr		# Cosmic ray list pointer
+int	col		# Cofluxn
+int	line		# Line
+real	flux		# Luminosity
+real	ratio		# Ratio
+real	wt		# Weight
+real	replace		# Sky value
+int	flag		# Flag value
+
+int	ncr
+errchk	realloc
+
+begin
+	if (CR_NCR(cr) == CR_NALLOC(cr)) {
+	    CR_NALLOC(cr) = CR_NALLOC(cr) + CR_ALLOC
+	    call realloc (CR_COL(cr), CR_NALLOC(cr), TY_REAL)
+	    call realloc (CR_LINE(cr), CR_NALLOC(cr), TY_REAL)
+	    call realloc (CR_FLUX(cr), CR_NALLOC(cr), TY_REAL)
+	    call realloc (CR_RATIO(cr), CR_NALLOC(cr), TY_REAL)
+	    call realloc (CR_WT(cr), CR_NALLOC(cr), TY_REAL)
+	    call realloc (CR_REPLACE(cr), CR_NALLOC(cr), TY_REAL)
+	    call realloc (CR_FLAG(cr), CR_NALLOC(cr), TY_INT)
+	}
+
+	ncr = CR_NCR(cr)
+	CR_NCR(cr) = ncr + 1
+	Memr[CR_COL(cr)+ncr] = col
+	Memr[CR_LINE(cr)+ncr] = line
+	Memr[CR_FLUX(cr)+ncr] = flux
+	Memr[CR_RATIO(cr)+ncr] = ratio
+	Memr[CR_WT(cr)+ncr] = wt
+	Memr[CR_REPLACE(cr)+ncr] = replace
+	Memi[CR_FLAG(cr)+ncr] = flag
+end
+
+
+# CR_TRAIN -- Set flux ratio threshold from a training set.
+
+procedure cr_train (cr, gp, gt, im, fluxratio, fname)
+
+pointer	cr		#I Cosmic ray list
+pointer	gp		#I GIO pointer
+pointer	gt		#I GTOOLS pointer
+pointer	im		#I IMIO pointer
+real	fluxratio	#O Flux ratio threshold
+char	fname[ARB]	#I Save file name
+
+char	cmd[10]
+bool	gflag
+real	x, y, y1, y2, w, r, rmin
+int	i, j, n, f, ncr, wcs, key, fd, clgcur(), open(), errcode()
+pointer	col, line, ratio, flux, wt, flag
+
+begin
+	# Open save file
+	iferr (fd = open (fname, APPEND, TEXT_FILE)) {
+	    if (errcode() != SYS_FNOFNAME)
+		call erract (EA_WARN)
+	    fd = 0
+	}
+
+	ncr = CR_NCR(cr)
+	col = CR_COL(cr) - 1
+	line = CR_LINE(cr) - 1
+	flux = CR_FLUX(cr) - 1
+	ratio = CR_RATIO(cr) - 1
+	wt = CR_WT(cr) - 1
+	flag = CR_FLAG(cr) - 1
+
+	gflag = false
+	n = 0
+	while (clgcur ("objects", x, y, wcs, key, cmd, 10) != EOF) {
+	    switch (key) {
+	    case '?':
+		call gpagefile (gp, HELP, PROMPT)
+		next
+	    case 'q':
+		break
+	    case 's':
+		w = 1
+		f = ALWAYSNO
+	    case 'c':
+		w = -1
+		f = ALWAYSYES
+	    case 'g':
+		if (gflag)
+		    call cr_examine (cr, gp, gt, im, fluxratio, 'z')
+		else {
+		    if (n > 1)
+			call cr_findthresh (cr, fluxratio)
+		    call cr_flags (cr, fluxratio)
+		    call cr_examine (cr, gp, gt, im, fluxratio, 'r')
+		    gflag = true
+		}
+		next
+	    default:
+		next
+	    }
+
+	    y1 = y - CR_RMAX
+	    y2 = y + CR_RMAX
+	    for (i=10; i<ncr && y1>Memr[line+i]; i=i+10)
+		;
+	    j = i - 9
+	    rmin = (Memr[col+j] - x) ** 2 + (Memr[line+j] - y) ** 2
+	    for (i=j+1; i<ncr && y2>Memr[line+i]; i=i+1) {
+		r = (Memr[col+i] - x) ** 2 + (Memr[line+i] - y) ** 2
+		if (r < rmin) {
+		    rmin = r
+		    j = i
+		}
+	    }
+	    if (sqrt (rmin) > CR_RMAX)
+		next
+
+	    Memr[wt+j] = w
+	    Memi[flag+j] = f
+	    n = n + 1
+
+	    if (gflag) {
+		if (n > 1) {
+		    call cr_findthresh (cr, r)
+		    call cr_update (gp, r, cr, fluxratio, 0)
+		}
+		call gmark (gp, Memr[flux+j], Memr[ratio+j], GM_BOX, 2., 2.)
+	    }
+	    if (fd > 0) {
+		call fprintf (fd, "%g %g %d %c\n")
+		    call pargr (x)
+		    call pargr (y)
+		    call pargi (wcs)
+		    call pargi (key)
+	    }
+	}
+
+	if (fd > 0)
+	    call close (fd)
+end
+
+
+# CR_FINDTHRESH -- Find flux ratio.
+
+procedure cr_findthresh (cr, fluxratio)
+
+pointer	cr		#I Cosmic ray list
+real	fluxratio	#O Flux ratio threshold
+
+real	w, r, rmin, cr_weight()
+int	i, ncr
+pointer	ratio, wt
+
+begin
+	ncr = CR_NCR(cr)
+	ratio = CR_RATIO(cr) - 1
+	wt = CR_WT(cr) - 1
+
+	fluxratio = Memr[ratio+1]
+	rmin = cr_weight (fluxratio, Memr[ratio+1], Memr[wt+1], ncr)
+	do i = 2, ncr {
+	    if (Memr[wt+i] == 0.)
+		next
+	    r = Memr[ratio+i]
+	    w = cr_weight (r, Memr[ratio+1], Memr[wt+1], ncr)
+	    if (w <= rmin) {
+		if (w == rmin)
+		    fluxratio = min (fluxratio, r)
+		else {
+		    rmin = w
+		    fluxratio = r
+		}
+	    }
+	}
+end
+
+
+# CR_WEIGHT -- Compute the training weight at a particular flux ratio.
+
+real procedure cr_weight (fluxratio, ratio, wts, ncr)
+
+real	fluxratio		#I Flux ratio
+real	ratio[ARB]		#I Ratio Values
+real	wts[ARB]		#I Weights
+int	ncr			#I Number of ratio values
+real	wt			#O Sum of weights
+
+int	i
+
+begin
+	wt = 0.
+	do i = 1, ncr {
+	    if (ratio[i] > fluxratio) {
+		if (wts[i] < 0.)
+		    wt = wt - wts[i]
+	    } else {
+		if (wts[i] > 0.)
+		    wt = wt + wts[i]
+	    }
+	}
+	return (wt)
+end
+
+
+# CR_FLAGS -- Set cosmic ray reject flags.
+
+procedure cr_flags (cr, fluxratio)
+
+pointer	cr			# Cosmic ray candidate list
+real	fluxratio		# Rejection limits
+
+int	i, ncr
+pointer	ratio, flag
+
+begin
+	ncr = CR_NCR(cr)
+	ratio = CR_RATIO(cr) - 1
+	flag = CR_FLAG(cr) - 1
+
+	do i = 1, ncr {
+	    if ((Memi[flag+i] == ALWAYSYES) || (Memi[flag+i] == ALWAYSNO))
+		next
+	    if (Memr[ratio+i] > fluxratio)
+		Memi[flag+i] = NO
+	    else
+		Memi[flag+i] = YES
+	}
+end
+
+
+# CR_BADPIX -- Store cosmic rays in bad pixel list.
+# This is currently a temporary measure until a real bad pixel list is
+# implemented.
+
+procedure cr_badpix (cr, fname)
+
+pointer	cr		# Cosmic ray list
+char	fname[ARB]	# Bad pixel file name
+
+int	i, ncr, c, l, f, fd, open(), errcode()
+pointer	col, line, ratio, flux, flag
+errchk	open
+
+begin
+	# Open bad pixel file
+	iferr (fd = open (fname, APPEND, TEXT_FILE)) {
+	    if (errcode() != SYS_FNOFNAME)
+		call erract (EA_WARN)
+	    return
+	}
+
+	ncr = CR_NCR(cr)
+	col = CR_COL(cr) - 1
+	line = CR_LINE(cr) - 1
+	flux = CR_FLUX(cr) - 1
+	ratio = CR_RATIO(cr) - 1
+	flag = CR_FLAG(cr) - 1
+
+	do i = 1, ncr {
+	    f = Memi[flag+i]
+	    if ((f == NO) || (f == ALWAYSNO))
+		next
+
+	    c = Memr[col+i]
+	    l = Memr[line+i]
+	    call fprintf (fd, "%d %d\n")
+	        call pargi (c)
+	        call pargi (l)
+	}
+	call close (fd)
+end
+
+
+# CR_CRMASK -- Store cosmic rays in cosmic ray mask.
+
+procedure cr_crmask (cr, fname, ref)
+
+pointer	cr		# Cosmic ray list
+char	fname[ARB]	# Cosmic ray mask
+pointer	ref		# Refrence image
+
+int	i, axlen[7], depth, ncr, f, rl[3,2]
+long	v[2]
+pointer	sp, title, pm, col, line, flag, pm_newmask()
+errchk	pm_loadf, pm_newmask, pm_savef
+
+begin
+	call smark (sp)
+	call salloc (title, SZ_LINE, TY_CHAR)
+
+	pm = pm_newmask (ref, 1)
+	iferr (call pm_loadf (pm, fname, Memc[title], SZ_LINE))
+	    ;
+	call pm_gsize (pm, i, axlen, depth)
+
+	ncr = CR_NCR(cr)
+	col = CR_COL(cr) - 1
+	line = CR_LINE(cr) - 1
+	flag = CR_FLAG(cr) - 1
+
+	RL_LEN(rl) = 2
+	RL_AXLEN(rl) = axlen[1]
+	RL_N(rl,2) = 1
+	RL_V(rl,2) = 1
+	RL_X(rl,2) = 1
+
+	do i = 1, ncr {
+	    f = Memi[flag+i]
+	    if ((f == NO) || (f == ALWAYSNO))
+		next
+
+	    v[1] = Memr[col+i]
+	    v[2] = Memr[line+i]
+	    call pmplri (pm, v, rl, depth, 1, PIX_SRC)
+	}
+
+	call strcpy ("Cosmic ray mask from COSMICRAYS", Memc[title], SZ_LINE)
+	call pm_savef (pm, fname, Memc[title], PM_UPDATE)
+	call pm_close (pm)
+	call sfree (sp)
+end
+
+
+# CR_REPLACE -- Replace cosmic rays in image with replacement values.
+
+procedure cr_replace (cr, offset, im, nreplaced)
+
+pointer	cr		# Cosmic ray list
+int	offset		# Offset in list
+pointer	im		# IMIO pointer of output image
+int	nreplaced	# Number replaced (for log)
+
+int	i, ncr, c, l, f
+real	r
+pointer	col, line, replace, flag, imps2r()
+
+begin
+	ncr = CR_NCR(cr)
+	if (ncr <= offset)
+	    return
+
+	col = CR_COL(cr) - 1
+	line = CR_LINE(cr) - 1
+	replace = CR_REPLACE(cr) - 1
+	flag = CR_FLAG(cr) - 1
+
+	do i = offset+1, ncr {
+	    f = Memi[flag+i]
+	    if ((f == NO) || (f == ALWAYSNO))
+		next
+
+	    c = Memr[col+i]
+	    l = Memr[line+i]
+	    r = Memr[replace+i]
+	    Memr[imps2r (im, c, c, l, l)] = r
+	    nreplaced = nreplaced + 1
+	}
+end
diff --git a/noao/imred/crutil/src/crmedian.par b/noao/imred/crutil/src/crmedian.par
new file mode 100644
index 00000000..2baccb05
--- /dev/null
+++ b/noao/imred/crutil/src/crmedian.par
@@ -0,0 +1,15 @@
+input,f,a,,,,Input image
+output,f,a,,,,Output image
+crmask,f,h,,,,Output cosmic ray mask
+median,f,h,,,,Output median image
+sigma,f,h,,,,Output sigma image
+residual,f,h,,,,Output residual image
+var0,r,h,0.,0.,,Variance coefficient for DN^0 term
+var1,r,h,0.,0.,,Variance coefficient for DN^1 term
+var2,r,h,0.,0.,,Variance coefficient for DN^2 term
+lsigma,r,h,10.,,,Low clipping sigma factor
+hsigma,r,h,3.,,,High clipping sigma factor
+ncmed,i,h,5,1,,Column box size for median level calculation
+nlmed,i,h,5,1,,Line box size for median level calculation
+ncsig,i,h,25,10,,Column box size for sigma calculation
+nlsig,i,h,25,10,,Line box size for sigma calculation
diff --git a/noao/imred/crutil/src/crnebula.cl b/noao/imred/crutil/src/crnebula.cl
new file mode 100644
index 00000000..c5f214f9
--- /dev/null
+++ b/noao/imred/crutil/src/crnebula.cl
@@ -0,0 +1,116 @@
+# CRNEBULA -- Cosmic ray cleaning for images with fine nebular structure.
+
+procedure crnebula (input, output)
+
+file    input           {prompt="Input image"}
+file	output		{prompt="Output image"}
+file    crmask          {prompt="Cosmic ray mask"}
+file	residual	{prompt="Residual median image"}
+file	rmedresid	{prompt="Residual ring median image"}
+real    var0 = 1.       {prompt="Variance coefficient for DN^0 term", min=0.}
+real	var1 = 0.	{prompt="Variance coefficient for DN^1 term", min=0.}
+real	var2 = 0.	{prompt="Variance coefficient for DN^2 term", min=0.}
+real    sigmed = 3      {prompt="Sigma clip factor for residual median"}
+real    sigrmed = 3     {prompt="Sigma clip factor for residual ring median"}
+
+int     mbox = 5        {prompt="Median box size"}
+real    rin = 1.5       {prompt="Inner radius for ring median"}
+real    rout = 6        {prompt="Outer radius for ring median"}
+bool    verbose = no    {prompt="Verbose"}
+
+begin
+        file    in, out, med, sig, resid, rmed
+        struct  expr
+
+        # Query once for query parameters.
+        in = input
+        out = output
+
+        # Check on output images.
+	if (out != "" && imaccess (out))
+            error (1, "Output image already exists ("//out//")")
+        if (crmask != "" && imaccess (crmask))
+            error (1, "Output mask already exists ("//crmask//")")
+	if (residual != "" && imaccess (residual))
+            error (1, "Output residual image already exists ("//residual//")")
+	if (rmedresid != "" && imaccess (rmedresid))
+            error (1,
+	       "Output ring median difference already exists ("//rmedresid//")")
+
+        # Create median results.
+	med = mktemp ("cr")
+	sig = mktemp ("cr")
+	resid = residual
+	if (resid == "")
+	    resid = mktemp ("cr")
+	if (verbose)
+	    printf ("Creating CRMEDIAN results\n")
+	crmedian (in, "", crmask="", median=med, sigma=sig, residual=resid,
+	    var0=var0, var1=var1, var2=var2, lsigma=100., hsigma=sigmed,
+	    ncmed=mbox, nlmed=mbox, ncsig=25, nlsig=25)
+
+        # Create ring median filtered image.
+        rmed = mktemp ("cr")
+	rmedian (in, rmed, rin, rout, ratio=1., theta=0., zloreject=INDEF,
+	    zhireject=INDEF, boundary="wrap", constant=0., verbose=verbose)
+
+	# Create output images.
+	if (rmedresid != "") {
+	    printf ("(a-b)/c\n") | scan (expr)
+	    imexpr (expr, rmedresid, med, rmed, sig, dims="auto",
+		intype="auto", outtype="real", refim="auto", bwidth=0,
+		btype="nearest", bpixval=0., rangecheck=yes, verbose=no,
+		exprdb="none")
+	    imdelete (rmed, verify-)
+	    imdelete (sig, verify-)
+	    if (out != "") {
+		if (verbose)
+		    printf ("Create output image %s\n", out)
+		printf ("((a<%.3g)||(abs(b)>%.3g)) ? c : d\n",
+		    sigmed, sigrmed) | scan (expr)
+		imexpr (expr, out, resid, rmedresid, in, med, dims="auto",
+		    intype="auto", outtype="auto", refim="auto", bwidth=0,
+		    btype="nearest", bpixval=0., rangecheck=yes, verbose=no,
+		    exprdb="none")
+	    }
+	    if (crmask != "") {
+		if (verbose)
+		    printf ("Create cosmic ray mask %s\n", crmask)
+		printf ("((a<%.3g)||(abs(b)>%.3g)) ? 0 : 1\n",
+		    sigmed, sigrmed) | scan (expr)
+		set imtype = "pl"
+		imexpr (expr, crmask, resid, rmedresid, dims="auto",
+		    intype="auto", outtype="short", refim="auto", bwidth=0,
+		    btype="nearest", bpixval=0., rangecheck=yes, verbose=no,
+		    exprdb="none")
+	    }
+	    imdelete (med, verify-)
+	} else {
+	    if (out != "") {
+		if (verbose)
+		    printf ("Create output image %s\n", out)
+		printf ("((a<%.3g)||(abs((b-c)/d)>%.3g)) ? e : b\n",
+		    sigmed, sigrmed) | scan (expr)
+		imexpr (expr, out, resid, med, rmed, sig, in, dims="auto",
+		    intype="auto", outtype="auto", refim="auto", bwidth=0,
+		    btype="nearest", bpixval=0., rangecheck=yes, verbose=no,
+		    exprdb="none")
+	    }
+	    if (crmask != "") {
+		if (verbose)
+		    printf ("Create cosmic ray mask %s\n", crmask)
+		printf ("((a<%.3g)||(abs((b-c)/d)>%.3g)) ? 0 : 1\n",
+		    sigmed, sigrmed) | scan (expr)
+		set imtype = "pl"
+		imexpr (expr, crmask, resid, med, rmed, sig, dims="auto",
+		    intype="auto", outtype="short", refim="auto", bwidth=0,
+		    btype="nearest", bpixval=0., rangecheck=yes, verbose=no,
+		    exprdb="none")
+	    }
+	    imdelete (med, verify-)
+	    imdelete (sig, verify-)
+	    imdelete (rmed, verify-)
+	}
+	if (residual == "")
+	    imdelete (resid, verify-)
+end
diff --git a/noao/imred/crutil/src/crsurface.x b/noao/imred/crutil/src/crsurface.x
new file mode 100644
index 00000000..32645ff4
--- /dev/null
+++ b/noao/imred/crutil/src/crsurface.x
@@ -0,0 +1,46 @@
+define	DUMMY	6
+
+# CR_SURFACE -- Draw a perspective view of a surface.  The altitude
+# and azimuth of the viewing angle are variable.
+
+procedure cr_surface(gp, data, ncols, nlines, angh, angv)
+
+pointer	gp			# GIO pointer
+real	data[ncols,nlines]	# Surface data to be plotted
+int	ncols, nlines		# Dimensions of surface
+real	angh, angv		# Orientation of surface (degrees)
+
+int	wkid
+pointer	sp, work
+
+int	first
+real	vpx1, vpx2, vpy1, vpy2
+common  /frstfg/ first
+common  /noaovp/ vpx1, vpx2, vpy1, vpy2
+
+begin
+	call smark (sp)
+	call salloc (work, 2 * (2 * ncols * nlines + ncols + nlines), TY_REAL)
+
+	# Initialize surface common blocks
+	first = 1
+	call srfabd()
+
+	# Define viewport.
+	call ggview (gp, vpx1, vpx2, vpy1, vpy2) 
+
+	# Link GKS to GIO
+	wkid = 1
+	call gopks (STDERR)
+	call gopwk (wkid, DUMMY, gp)
+	call gacwk (wkid)
+
+	call ezsrfc (data, ncols, nlines, angh, angv, Memr[work])
+
+	call gdawk (wkid)
+	# We don't want to close the GIO pointer.
+	#call gclwk (wkid)
+	call gclks ()
+
+	call sfree (sp)
+end
diff --git a/noao/imred/crutil/src/mkpkg b/noao/imred/crutil/src/mkpkg
new file mode 100644
index 00000000..1279b656
--- /dev/null
+++ b/noao/imred/crutil/src/mkpkg
@@ -0,0 +1,38 @@
+# COSMIC RAY CLEANING PACKAGE
+
+$call	relink
+$exit
+
+update:
+        $call   relink
+        $call   install
+        ;
+
+relink:
+        $update libpkg.a
+        $call   crutil
+        ;
+
+install:
+        $move   xx_crutil.e noaobin$x_crutil.e
+        ;
+
+crutil:
+        $omake  x_crutil.x
+        $link   x_crutil.o libpkg.a -lxtools -lcurfit -lgsurfit -lncar -lgks\
+                -o xx_crutil.e
+        ;
+
+libpkg.a:
+	crexamine.x	crlist.h <error.h> <gset.h> <mach.h> <pkg/gtools.h>\
+			<imhdr.h>
+	crfind.x	<math/gsurfit.h>
+	crlist.x	crlist.h <error.h> <gset.h> <pmset.h>
+	crsurface.x	
+	t_cosmicrays.x	crlist.h <error.h> <math/gsurfit.h> <imhdr.h> <gset.h>\
+			<pkg/gtools.h> <imset.h>
+	t_craverage.x	<imhdr.h> <error.h> <mach.h>
+	t_crgrow.x	<error.h> <imhdr.h>
+	t_crmedian.x	<imhdr.h> <mach.h>
+	xtmaskname.x
+	;
diff --git a/noao/imred/crutil/src/t_cosmicrays.x b/noao/imred/crutil/src/t_cosmicrays.x
new file mode 100644
index 00000000..fa68a39d
--- /dev/null
+++ b/noao/imred/crutil/src/t_cosmicrays.x
@@ -0,0 +1,329 @@
+include	<error.h>
+include <imhdr.h>
+include <imset.h>
+include	<math/gsurfit.h>
+include	<gset.h>
+include	<pkg/gtools.h>
+include	"crlist.h"
+
+# T_COSMICRAYS -- Detect and remove cosmic rays in images.
+# A list of images is examined for cosmic rays which are then replaced
+# by values from neighboring pixels.  The output image may be the same
+# as the input image.  This is the top level procedure which manages
+# the input and output image data.  The actual algorithm for detecting
+# cosmic rays is in CR_FIND.
+
+procedure t_cosmicrays ()
+
+int	list1			# List of input images to be cleaned
+int	list2			# List of output images
+int	list3			# List of output bad pixel files
+real	threshold		# Detection threshold
+real	fluxratio		# Luminosity boundary for stars
+int	npasses			# Number of cleaning passes
+int	szwin			# Size of detection window
+bool	train			# Use training objects?
+pointer	savefile		# Save file for training objects
+bool	interactive		# Examine cosmic ray parameters?
+char	ans			# Answer to interactive query
+
+int	nc, nl, c, c1, c2, l, l1, l2, szhwin, szwin2
+int	i, j, k, m, ncr, ncrlast, nreplaced, flag
+pointer	sp, input, output, badpix, str, gp, gt, im, in, out
+pointer	x, y, z, w, sf1, sf2, cr, data, ptr
+
+bool	clgetb(), streq(), strne()
+char	clgetc()
+int	imtopenp(), imtlen(), imtgetim(), clpopnu(), clgfil(), clgeti()
+real	clgetr()
+pointer	immap(), impl2r(), imgs2r(), gopen(), gt_init()
+errchk	immap, impl2r, imgs2r
+errchk	cr_find, cr_examine, cr_replace, cr_plot, cr_crmask
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (badpix, SZ_FNAME, TY_CHAR)
+	call salloc (savefile, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Get the task parameters.  Check that the number of output images
+	# is either zero, in which case the cosmic rays will be removed
+	# in place, or equal to the number of input images.
+
+	list1 = imtopenp ("input")
+	list2 = imtopenp ("output")
+	i = imtlen (list1)
+	j = imtlen (list2)
+	if (j > 0 && j != i)
+	    call error (0, "Input and output image lists do not match")
+
+	list3 = clpopnu ("crmasks")
+	threshold = clgetr ("threshold")
+	fluxratio = clgetr ("fluxratio")
+	npasses = clgeti ("npasses")
+	szwin = clgeti ("window")
+	train = clgetb ("train")
+	call clgstr ("savefile", Memc[savefile], SZ_FNAME)
+	interactive = clgetb ("interactive")
+	call clpstr ("answer", "yes")
+	ans = 'y'
+
+	# Set up the graphics.
+	call clgstr ("graphics", Memc[str], SZ_LINE)
+	if (interactive) {
+	    gp = gopen (Memc[str], NEW_FILE+AW_DEFER, STDGRAPH)
+	    gt = gt_init()
+	    call gt_sets (gt, GTTYPE, "mark")
+	    call gt_sets (gt, GTXTRAN, "log")
+	    call gt_setr (gt, GTXMIN, 10.)
+	    call gt_setr (gt, GTYMIN, 0.)
+	    call gt_sets (gt, GTTITLE, "Parameters of cosmic rays candidates")
+	    call gt_sets (gt, GTXLABEL, "Flux")
+	    call gt_sets (gt, GTYLABEL, "Flux Ratio")
+	}
+
+	# Set up surface fitting.  The background points are placed together
+	# at the beginning of the arrays.  There are two surface pointers,
+	# one for using the fast refit if there are no points excluded and
+	# one for doing a full fit with points excluded.
+
+	szhwin = szwin / 2
+	szwin2 = szwin * szwin
+	call salloc (data, szwin, TY_INT)
+	call salloc (x, szwin2, TY_REAL)
+	call salloc (y, szwin2, TY_REAL)
+	call salloc (z, szwin2, TY_REAL)
+	call salloc (w, szwin2, TY_REAL)
+
+	k = 0
+	do i = 1, szwin {
+	    Memr[x+k] = i
+	    Memr[y+k] = 1
+	    k = k + 1
+	}
+	do i = 2, szwin {
+	    Memr[x+k] = szwin
+	    Memr[y+k] = i
+	    k = k + 1
+	}
+	do i = szwin-1, 1, -1 {
+	    Memr[x+k] = i
+	    Memr[y+k] = szwin
+	    k = k + 1
+	}
+	do i = szwin-1, 2, -1 {
+	    Memr[x+k] = 1
+	    Memr[y+k] = i
+	    k = k + 1
+	}
+	do i = 2, szwin-1 {
+	    do j = 2, szwin-1 {
+	        Memr[x+k] = j
+	        Memr[y+k] = i
+	        k = k + 1
+	    }
+	}
+	call aclrr (Memr[z], szwin2)
+	call amovkr (1., Memr[w], 4*szwin-4)
+	call gsinit (sf1, GS_POLYNOMIAL, 2, 2, NO, 1., real(szwin),
+	    1., real(szwin))
+	call gsinit (sf2, GS_POLYNOMIAL, 2, 2, NO, 1., real(szwin),
+	    1., real(szwin))
+	call gsfit (sf1, Memr[x], Memr[y], Memr[z], Memr[w], 4*szwin-4,
+	    WTS_USER, j)
+
+	# Process each input image.  Either work in place or create a
+	# new output image.  If an error mapping the images occurs
+	# issue a warning and go on to the next input image.
+
+	while (imtgetim (list1, Memc[input], SZ_FNAME) != EOF) {
+	    if (imtgetim (list2, Memc[output], SZ_FNAME) == EOF)
+		call strcpy (Memc[input], Memc[output], SZ_FNAME)
+	    if (clgfil (list3, Memc[badpix], SZ_FNAME) == EOF)
+		Memc[badpix] = EOS
+
+	    iferr {
+		in = NULL
+		out = NULL
+		cr = NULL
+
+		# Map the input image.  If the output image is
+		# the same as the input image work in place.
+		# Initialize IMIO to use a scrolling buffer of lines.
+
+	        if (streq (Memc[input], Memc[output])) {
+	            im = immap (Memc[input], READ_WRITE, 0)
+		} else
+		    im = immap (Memc[input], READ_ONLY, 0)
+		in = im
+
+	        nc = IM_LEN(in,1)
+	        nl = IM_LEN(in,2)
+	        if ((nl < szwin) || (nc < szwin))
+	            call error (0, "Image size is too small")
+	        call imseti (in, IM_NBUFS, szwin)
+	        call imseti (in, IM_TYBNDRY, BT_NEAREST)
+	        call imseti (in, IM_NBNDRYPIX, szhwin)
+
+		# Open the output image if needed.
+	        if (strne (Memc[input], Memc[output]))
+		    im = immap (Memc[output], NEW_COPY, in)
+		out = im
+
+		# Open a cosmic ray list structure.
+	        call cr_open (cr)
+		ncrlast = 0
+		nreplaced = 0
+
+		# Now proceed through the image line by line, scrolling
+		# the line buffers at each step.  If creating a new image
+		# also write out each line as it is read.  A procedure is
+		# called to find the cosmic ray candidates in the line
+		# and add them to the list maintained by CRLIST.
+		# Note that cosmic rays are not replaced at this point
+		# in order to allow the user to modify the criteria for
+		# a cosmic ray and review the results.
+
+	        c1 = 1-szhwin
+	        c2 = nc+szhwin
+		do i = 1, szwin-1 
+		    Memi[data+i] =
+			imgs2r (in, c1, c2, i-szhwin, i-szhwin)
+
+	        do l = 1, nl {
+		    do i = 1, szwin-1
+			Memi[data+i-1] = Memi[data+i]
+		    Memi[data+szwin-1] =
+			imgs2r (in, c1, c2, l+szhwin, l+szhwin)
+		    if (out != in)
+		        call amovr (Memr[Memi[data+szhwin]+szhwin],
+			    Memr[impl2r(out,l)], nc)
+
+	            call cr_find (cr, threshold, Memi[data],
+		        c2-c1+1, szwin, c1, l,
+		        sf1, sf2, Memr[x], Memr[y], Memr[z], Memr[w])
+	        }
+		if (interactive && train) {
+		    call cr_train (cr, gp, gt, in, fluxratio, Memc[savefile])
+		    train = false
+		}
+	        call cr_flags (cr, fluxratio)
+
+		# If desired examine the cosmic ray list interactively.
+	        if (interactive && ans != 'N') {
+		    if (ans != 'Y') {
+		        call eprintf ("%s - ")
+		            call pargstr (Memc[input])
+		        call flush (STDERR)
+		        ans = clgetc ("answer")
+		    }
+		    if ((ans == 'Y') || (ans == 'y'))
+	                call cr_examine (cr, gp, gt, in, fluxratio, 'r')
+	        }
+
+		# Now replace the selected cosmic rays in the output image.
+
+		call imflush (out)
+		call imseti (out, IM_ADVICE, RANDOM)
+	        call cr_replace (cr, ncrlast, out, nreplaced)
+
+		# Do additional passes through the data.  We work in place
+		# in the output image.  Note that we only have to look in
+		# the vicinity of replaced cosmic rays for secondary
+		# events since we've already looked at every pixel once.
+		# Instead of scrolling through the image we will extract
+		# subrasters around each replaced cosmic ray.  However,
+		# we use pointers into the subraster to maintain the same
+		# format expected by CR_FIND.
+
+	        if (npasses > 1) {
+	            if (out != in)
+	                call imunmap (out)
+	            call imunmap (in)
+		    im = immap (Memc[output], READ_WRITE, 0)
+		    in = im
+		    out = im
+	            call imseti (in, IM_TYBNDRY, BT_NEAREST)
+	            call imseti (in, IM_NBNDRYPIX, szhwin)
+
+	            for (i=2; i<=npasses; i=i+1) {
+			# Loop through each cosmic ray in the previous pass.
+		        ncr = CR_NCR(cr)
+		        do j = ncrlast+1, ncr {
+			    flag = Memi[CR_FLAG(cr)+j-1]
+			    if (flag==NO || flag==ALWAYSNO)
+			        next
+			    c = Memr[CR_COL(cr)+j-1]
+			    l = Memr[CR_LINE(cr)+j-1]
+			    c1 = max (1-szhwin, c - (szwin-1))
+			    c2 = min (nc+szhwin, c + (szwin-1))
+			    k = c2 - c1 + 1
+			    l1 = max (1-szhwin, l - (szwin-1))
+			    l2 = min (nl+szhwin, l + (szwin-1))
+
+			    # Set the line pointers off an image section
+			    # centered on a previously replaced cosmic ray.
+
+			    ptr = imgs2r (in, c1, c2, l1, l2) - k
+
+			    l1 = max (1, l - szhwin)
+			    l2 = min (nl, l + szhwin)
+			    do l = l1, l2 {
+				do m = 1, szwin
+				    Memi[data+m-1] = ptr + m * k
+				ptr = ptr + k
+
+	                        call cr_find ( cr, threshold, Memi[data],
+				    k, szwin, c1, l, sf1, sf2,
+				    Memr[x], Memr[y], Memr[z], Memr[w])
+			    }
+	                }
+		        call cr_flags (cr, fluxratio)
+
+			# Replace any new cosmic rays found.
+	                call cr_replace (cr, ncr, in, nreplaced)
+			ncrlast = ncr
+		    }
+	        }
+
+		# Output header log, log, plot, and bad pixels.
+		call sprintf (Memc[str], SZ_LINE,
+		    "Threshold=%5.1f, fluxratio=%6.2f, removed=%d")
+		    call pargr (threshold)
+		    call pargr (fluxratio)
+		    call pargi (nreplaced)
+		call imastr (out, "crcor", Memc[str])
+
+		call cr_plot (cr, in, fluxratio)
+		call cr_crmask (cr, Memc[badpix], in)
+
+	        call cr_close (cr)
+	        if (out != in)
+	            call imunmap (out)
+	        call imunmap (in)
+	    } then {
+		# In case of error clean up and go on to the next image.
+		if (in != NULL) {
+		    if (out != NULL && out != in)
+			call imunmap (out)
+		    call imunmap (in)
+		}
+		if (cr != NULL)
+		    call cr_close (cr)
+		call erract (EA_WARN)
+	    }
+	}
+
+	if (interactive) {
+	    call gt_free (gt)
+	    call gclose (gp)
+	}
+	call imtclose (list1)
+	call imtclose (list2)
+	call clpcls (list3)
+	call gsfree (sf1)
+	call gsfree (sf2)
+	call sfree (sp)
+end
diff --git a/noao/imred/crutil/src/t_craverage.x b/noao/imred/crutil/src/t_craverage.x
new file mode 100644
index 00000000..f7b82113
--- /dev/null
+++ b/noao/imred/crutil/src/t_craverage.x
@@ -0,0 +1,847 @@
+include	<error.h>
+include	<imhdr.h>
+include	<mach.h>
+
+define	MAXBUF		500000	# Maximum pixel buffer
+
+define	PLSIG		15.87	# Low percentile
+define	PHSIG		84.13	# High percentile
+
+
+# T_CRAVERAGE -- Detect, fix, and flag cosmic rays.  Also detect objects.
+# Deviant pixels relative to a local average with the candidate pixel
+# excluded and sigma are detected and replaced by the average value
+# and/or written to a cosmic ray mask.  Average values above a the median
+# of a background annulus are detected as objects and cosmic rays are
+# excluded.  The object positions may be output in the mask.
+
+procedure t_craverage ()
+
+int	inlist			# Input image list
+int	outlist			# Output image list
+int	crlist			# Output mask list
+int	avglist			# Output average list
+int	siglist			# Output sigma list
+int	crval			# Output cosmic ray mask value
+int	objval			# Output object mask value
+int	navg			# Averaging box size
+int	nrej			# Number of high pixels to reject from average
+int	nbkg			# Background width
+int	nsig			# Sigma box size
+real	lobjsig, hobjsig	# Object threshold sigmas
+real	lcrsig, hcrsig		# CR threshold sigmas outside of object
+real	var0			# Variance coefficient for DN^0 term
+real	var1			# Variance coefficient for DN^1 term
+real	var2			# Variance coefficient for DN^2 term
+real	crgrw			# Cosmic ray grow radius
+real	objgrw			# Object grow radius
+
+int	i, nc, nl, nlstep, nbox, l1, l2, l3, l4, nl1, pmmode
+pointer	sp, input, output, crmask, crmask1, extname, average, sigma
+pointer	in, out, pm, aim, sim
+pointer	inbuf, pinbuf, outbuf, pbuf, abuf, sbuf
+
+real	clgetr()
+int	clgeti(), imtopenp(), imtgetim()
+pointer	immap(), imgs2s(), imgs2r(), imps2r(), imps2s()
+errchk	immap, imgs2s, imgs2r, imps2r, imps2s, craverage, crgrow, imgstr
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (crmask, SZ_FNAME, TY_CHAR)
+	call salloc (crmask1, SZ_FNAME, TY_CHAR)
+	call salloc (average, SZ_FNAME, TY_CHAR)
+	call salloc (sigma, SZ_FNAME, TY_CHAR)
+	call salloc (extname, SZ_FNAME, TY_CHAR)
+
+	# Get parameters.
+	inlist = imtopenp ("input")
+	outlist = imtopenp ("output")
+	crlist = imtopenp ("crmask")
+	avglist = imtopenp ("average")
+	siglist = imtopenp ("sigma")
+	crval = clgeti ("crval")
+	objval = clgeti ("objval")
+	navg = max (1, clgeti ("navg") / 2)
+	nrej = min (clgeti ("nrej"), navg-1)
+	nbkg = clgeti ("nbkg")
+	nsig = clgeti ("nsig")
+	lobjsig = clgetr ("lobjsig")
+	hobjsig = clgetr ("hobjsig")
+	lcrsig = clgetr ("lcrsig")
+	hcrsig = clgetr ("hcrsig")
+	nbox = 2 * (navg + nbkg) + 1
+	var0 = clgetr ("var0")
+	var1 = clgetr ("var1")
+	var2 = clgetr ("var2")
+	crgrw = clgetr ("crgrow")
+	objgrw = clgetr ("objgrow")
+
+	# Do the input images.
+	Memc[crmask1] = EOS
+	while (imtgetim (inlist, Memc[input], SZ_FNAME) != EOF) {
+	    if (imtgetim (outlist, Memc[output], SZ_FNAME) == EOF)
+		Memc[output] = EOS
+	    if (imtgetim (crlist, Memc[crmask], SZ_FNAME) == EOF)
+	        call strcpy (Memc[crmask1], Memc[crmask], SZ_FNAME)
+	    else if (Memc[crmask] == '!')
+	        call strcpy (Memc[crmask], Memc[crmask1], SZ_FNAME)
+	    if (imtgetim (avglist, Memc[average], SZ_FNAME) == EOF)
+		Memc[average] = EOS
+	    if (imtgetim (siglist, Memc[sigma], SZ_FNAME) == EOF)
+		Memc[sigma] = EOS
+
+	    # Map the input and output images.
+	    iferr {
+		in = NULL; out = NULL; pm = NULL; aim = NULL; sim = NULL
+		inbuf = NULL; pinbuf = NULL; outbuf = NULL; pbuf = NULL;
+		abuf = NULL; sbuf=NULL
+
+		in = immap (Memc[input], READ_ONLY, 0)
+		if (Memc[output] != EOS)
+		    out = immap (Memc[output], NEW_COPY, in)
+		if (Memc[crmask] != EOS) {
+		    if (Memc[crmask] == '!')
+			call imgstr (in, Memc[crmask+1], Memc[crmask], SZ_FNAME)
+		    pmmode = READ_WRITE
+		    iferr (call imgstr (in, "extname", Memc[extname], SZ_FNAME))
+			call strcpy ("pl", Memc[extname], SZ_FNAME)
+		    call xt_maskname (Memc[crmask], Memc[extname], pmmode,
+			Memc[crmask], SZ_FNAME)
+		    iferr (pm = immap (Memc[crmask], pmmode, 0)) {
+			pmmode = NEW_COPY
+			pm = immap (Memc[crmask], pmmode, in)
+		    }
+		}
+		if (Memc[average] != EOS)
+		    aim = immap (Memc[average], NEW_COPY, in)
+		if (Memc[sigma] != EOS)
+		    sim = immap (Memc[sigma], NEW_COPY, in)
+
+		# Go through the input in large blocks of lines.  If the
+		# block is smaller than the whole image overlap the blocks
+		# so the average only has boundaries at the ends of the image.
+		# However, the output is done in non-overlapping blocks with
+		# the pointers are adjusted so that addresses can be in the
+		# space of the input block.  CRAVERAGE does not address
+		# outside of the output data block.  Set the mask values
+		# based on the distances to the nearest good pixels.
+
+		nc = IM_LEN(in,1)
+		nl = IM_LEN(in,2)
+		nlstep = max (1, MAXBUF / nc - nbox)
+
+		do i = 1, nl, nlstep {
+		    l1 = i
+		    l2 = min (nl, i + nlstep - 1)
+		    l3 = max (1, l1 - nbox / 2)
+		    l4 = min (nl, l2 + nbox / 2)
+		    nl1 = l4 - l3 + 1
+		    inbuf = imgs2r (in, 1, nc, l3, l4)
+		    if (out != NULL)
+			outbuf = imps2r (out, 1, nc, l1, l2) - (l1 - l3) * nc
+		    if (pm != NULL) {
+			if (pmmode == READ_WRITE) {
+			    pinbuf = imgs2s (pm, 1, nc, l3, l4)
+			    pbuf = imps2s (pm, 1, nc, l1, l2)
+			    call amovs (Mems[pinbuf+(l1-l3)*nc],
+				Mems[pbuf], nc*(l2-l1+1))
+			    pbuf = pbuf - (l1 - l3) * nc
+			} else {
+			    pinbuf = NULL
+			    pbuf = imps2s (pm, 1, nc, l1, l2)
+			    call aclrs (Mems[pbuf], nc*(l2-l1+1))
+			    pbuf = pbuf - (l1 - l3) * nc
+			}
+		    }
+		    if (aim != NULL)
+			abuf = imps2r (aim, 1, nc, l1, l2) - (l1 - l3) * nc
+		    if (sim != NULL)
+			sbuf = imps2r (sim, 1, nc, l1, l2) - (l1 - l3) * nc
+		    if (pinbuf == NULL)
+			call craverage (inbuf, outbuf, pbuf, abuf, sbuf,
+			    nc, nl1, l1-l3+1, l2-l3+1, navg, nrej, nbkg,
+			    var0, var1, var2, nsig, lcrsig, hcrsig,
+			    lobjsig, hobjsig, crval, objval)
+		    else
+			call craverage1 (inbuf, pinbuf, outbuf, pbuf, abuf,
+			    sbuf, nc, nl1, l1-l3+1, l2-l3+1, navg, nrej, nbkg,
+			    var0, var1, var2, nsig, lcrsig, hcrsig,
+			    lobjsig, hobjsig, crval, objval)
+		}
+
+		# Grow regions if desired.  The routines are nops if the
+		# grow is zero.
+
+		if (pm != NULL) {
+		    if (pmmode != READ_WRITE) {
+			call imunmap (pm)
+			iferr (pm = immap (Memc[crmask], READ_WRITE, 0))
+			    call error (1, "Can't reopen mask for growing")
+		    }
+
+		    if (crval == objval)
+			call crgrow (pm, max (crgrw, objgrw), crval, crval)
+		    else {
+			call crgrow (pm, crgrw, crval, crval)
+			call crgrow (pm, objgrw, objval, objval)
+		    }
+		}
+	    } then
+		call erract (EA_WARN)
+
+	    if (sim != NULL)
+		call imunmap (sim)
+	    if (aim != NULL)
+		call imunmap (aim)
+	    if (pm != NULL)
+		call imunmap (pm)
+	    if (out != NULL)
+		call imunmap (out)
+	    call imunmap (in)
+	}
+
+	call imtclose (inlist)
+	call imtclose (outlist)
+	call imtclose (crlist)
+	call imtclose (avglist)
+	call imtclose (siglist)
+
+	call sfree (sp)
+end
+
+
+# CRAVERAGE -- Detect, replace, and flag cosmic rays.
+# A local background is computed using moving box averages to avoid
+# contaminating bad pixels.  If variance model is given then that is
+# used otherwise a local sigma is computed in blocks (it is not a moving box
+# for efficiency) by using a percentile point of the sorted pixel values to
+# estimate the width of the distribution uncontaminated by bad pixels).  Once
+# the background and sigma are known deviant pixels are found by using sigma
+# threshold factors.
+
+procedure craverage (in, out, pout, aout, sout, nc, nl, nl1, nl2,
+	navg, nrej, nbkg, var0, var1, var2, nsig, lcrsig, hcrsig,
+	lobjsig, hobjsig crval, objval)
+
+pointer	in			#I Input data
+pointer	out			#O Output data
+pointer	pout			#O Output mask (0=good, 1=bad)
+pointer	aout			#O Output averages
+pointer	sout			#O Output sigmas
+int	nc, nl			#I Number of columns and lines
+int	nl1, nl2		#I Lines to compute
+int	navg			#I Averaging box half-size
+int	nrej			#I Number of high pixels to reject from average
+int	nbkg			#I Median background width
+real	var0			#I Variance coefficient for DN^0 term
+real	var1			#I Variance coefficient for DN^1 term
+real	var2			#I Variance coefficient for DN^2 term
+int	nsig			#I Sigma box size
+real	lcrsig, hcrsig		#I Threshold sigmas outside of object
+real	lobjsig, hobjsig	#I Object threshold sigmas
+int	crval			#I CR mask value
+int	objval			#I Object mask value
+
+int	i, j, c, c1, c2, c3, c4, l, l1, l2, l3, l4, n1, n2
+int	navg2, nbkg2, nsig2, plsig, phsig
+real	data, avg, bkg, sigma, losig, hosig
+real	low, high, cravg(), amedr()
+pointer	stack, avgs, bkgs, sigs, work1, work2
+pointer	ptr1, ptr2, ip, op, pp, ap, sp
+
+begin
+	navg2 = (2 * navg + 1) ** 2
+	nbkg2 = (2 * (navg + nbkg) + 1) ** 2 - navg2
+	nsig2 = nsig * nsig
+
+	call smark (stack)
+	call salloc (avgs, nc, TY_REAL)
+	call salloc (bkgs, nc, TY_REAL)
+	call salloc (sigs, nc, TY_REAL)
+	call salloc (work1, navg2, TY_REAL)
+	call salloc (work2, max (nsig2, nbkg2), TY_REAL)
+
+	if (var0 != 0. && var1 == 0. && var2 ==0.)
+	    call amovkr (sqrt(var0), Memr[sigs], nc)
+
+	avgs = avgs - 1
+	sigs = sigs - 1
+	bkgs = bkgs - 1
+
+	plsig = nint (PLSIG*nsig2/100.-1)
+	phsig = nint (PHSIG*nsig2/100.-1)
+	losig = lobjsig / sqrt (real(navg2-1))
+	hosig = hobjsig / sqrt (real(navg2-1))
+
+	do l = nl1, nl2 {
+	    # Compute statistics.
+	    l1 = max (1, l-navg-nbkg)
+	    l2 = max (1, l-navg)
+	    l3 = min (nl, l+navg)
+	    l4 = min (nl, l+navg+nbkg)
+	    ap = aout + (l - 1) * nc
+	    do c = 1, nc {
+		c1 = max (1, c-navg-nbkg)
+		c2 = max (1, c-navg)
+		c3 = min (nc, c+navg)
+		c4 = min (nc, c+navg+nbkg)
+		ptr1 = work1
+		ptr2 = work2
+		n1 = 0
+		n2 = 0
+		do j = l1, l2-1 {
+		    ip = in + (j - 1) * nc + c1 - 1
+		    do i = c1, c4 {
+			Memr[ptr2] = Memr[ip]
+			n2 = n2 + 1
+			ptr2 = ptr2 + 1
+			ip = ip + 1
+		    }
+		}
+		do j = l2, l3 {
+		    ip = in + (j - 1) * nc + c1 - 1
+		    do i = c1, c2-1 {
+			Memr[ptr2] = Memr[ip]
+			n2 = n2 + 1
+			ptr2 = ptr2 + 1
+			ip = ip + 1
+		    }
+		    do i = c2, c3 {
+			if (j != l || i != c) {
+			    Memr[ptr1] = Memr[ip]
+			    n1 = n1 + 1
+			    ptr1 = ptr1 + 1
+			}
+			ip = ip + 1
+		    }
+		    do i = c3+1, c4 {
+			Memr[ptr2] = Memr[ip]
+			n2 = n2 + 1
+			ptr2 = ptr2 + 1
+			ip = ip + 1
+		    }
+		}
+		do j = l3+1, l4 {
+		    ip = in + (j - 1) * nc + c1 - 1
+		    do i = c1, c4 {
+			Memr[ptr2] = Memr[ip]
+			n2 = n2 + 1
+			ptr2 = ptr2 + 1
+			ip = ip + 1
+		    }
+		}
+		avg = cravg (Memr[work1], n1, nrej)
+		bkg = amedr (Memr[work2], n2)
+		Memr[bkgs+c] = bkg
+		Memr[avgs+c] = avg
+		if (aout != NULL) {
+		    Memr[ap] = avg - bkg
+		    ap = ap + 1
+		}
+	    }
+
+	    # Compute sigmas and output if desired.
+	    if (var0 != 0. || var1 != 0. || var2 != 0.) {
+		if (var1 != 0.) {
+		    if (var2 != 0.) {
+			do c = 1, nc {
+			    data = max (0., Memr[avgs+c])
+			    Memr[sigs+c] = sqrt (var0+var1*data+var2*data**2)
+			}
+		    } else {
+			do c = 1, nc {
+			    data = max (0., Memr[avgs+c])
+			    Memr[sigs+c] = sqrt (var0 + var1 * data)
+			}
+		    }
+		} else if (var2 != 0.) {
+		    do c = 1, nc {
+			data = max (0., Memr[avgs+c])
+			Memr[sigs+c] = sqrt (var0 + var2 * data**2)
+		    }
+		}
+	    } else {
+		# Compute sigmas from percentiles.  This is done in blocks.
+		if (mod (l-nl1, nsig) == 0 && l<nl-nsig+1) {
+		    do c = 1, nc-nsig+1, nsig {
+			ptr2 = work2
+			n2 = 0
+			do j = l, l+nsig-1 {
+			    ip = in + (j - 1) * nc + c - 1
+			    do i = 1, nsig {
+				Memr[ptr2] = Memr[ip]
+				n2 = n2 + 1
+				ptr2 = ptr2 + 1
+				ip = ip + 1
+			    }
+			}
+			call asrtr (Memr[work2], Memr[work2], n2)
+			sigma = (Memr[work2+phsig] - Memr[work2+plsig]) / 2.
+			call amovkr (sigma, Memr[sigs+c], nsig)
+		    }
+		    call amovkr (sigma, Memr[sigs+c], nc-c+1)
+		}
+	    }
+	    if (sout != NULL) {
+		sp = sout + (l - 1) * nc
+		do c = 1, nc {
+		    Memr[sp] = Memr[sigs+c]
+		    sp = sp + 1
+		}
+	    }
+
+	    # Detect, fix, and flag cosmic rays.
+	    if (pout == NULL && out == NULL)
+		;
+	    else if (pout == NULL) {
+		ip = in + (l - 1) * nc
+		op = out + (l - 1) * nc
+		do c = 1, nc {
+		    data = Memr[ip]
+		    avg = Memr[avgs+c]
+		    bkg = Memr[bkgs+c]
+		    sigma = Memr[sigs+c]
+		    low = bkg - losig * sigma
+		    high = bkg + hosig * sigma
+		    if (avg < low || avg > high) {
+			Memr[op] = data
+		    } else {
+			low = avg - lcrsig * sigma
+			high = avg + hcrsig * sigma
+			if (data < low || data > high)
+			    Memr[op] = avg
+			else
+			    Memr[op] = data
+		    }
+		    ip = ip + 1
+		    op = op + 1
+		}
+	    } else if (out == NULL) {
+		ip = in + (l - 1) * nc
+		pp = pout + (l - 1) * nc
+		do c = 1, nc {
+		    data = Memr[ip]
+		    avg = Memr[avgs+c]
+		    bkg = Memr[bkgs+c]
+		    sigma = Memr[sigs+c]
+		    low = bkg - losig * sigma
+		    high = bkg + hosig * sigma
+		    if (avg < low || avg > high)
+			Mems[pp] = objval
+		    else {
+			low = avg - lcrsig * sigma
+			high = avg + hcrsig * sigma
+			if (data < low || data > high)
+			    Mems[pp] = crval
+		    }
+		    ip = ip + 1
+		    pp = pp + 1
+		}
+	    } else {
+		ip = in + (l - 1) * nc
+		op = out + (l - 1) * nc
+		pp = pout + (l - 1) * nc
+		do c = 1, nc {
+		    data = Memr[ip]
+		    avg = Memr[avgs+c]
+		    bkg = Memr[bkgs+c]
+		    sigma = Memr[sigs+c]
+		    low = bkg - losig * sigma
+		    high = bkg + hosig * sigma
+		    if (avg < low || avg > high) {
+			Memr[op] = data
+			Mems[pp] = objval
+		    } else {
+			low = avg - lcrsig * sigma
+			high = avg + hcrsig * sigma
+			if (data < low || data > high) {
+			    Memr[op] = avg
+			    Mems[pp] = crval
+			} else
+			    Memr[op] = data
+		    }
+		    ip = ip + 1
+		    op = op + 1
+		    pp = pp + 1
+		}
+	    }
+	}
+
+	call sfree (stack)
+end
+
+
+# CRAVERAGE1 -- Detect, replace, and flag cosmic rays checking input mask.
+# A local background is computed using moving box averages to avoid
+# contaminating bad pixels.  If variance model is given then that is
+# used otherwise a local sigma is computed in blocks (it is not a moving box
+# for efficiency) by using a percentile point of the sorted pixel values to
+# estimate the width of the distribution uncontaminated by bad pixels).  Once
+# the background and sigma are known deviant pixels are found by using sigma
+# threshold factors.
+
+procedure craverage1 (in, pin, out, pout, aout, sout, nc, nl, nl1, nl2,
+	navg, nrej, nbkg, var0, var1, var2, nsig, lcrsig, hcrsig,
+	lobjsig, hobjsig crval, objval)
+
+pointer	in			#I Input data
+pointer	pin			#I Pixel mask data
+pointer	out			#O Output data
+pointer	pout			#O Output mask (0=good, 1=bad)
+pointer	aout			#O Output averages
+pointer	sout			#O Output sigmas
+int	nc, nl			#I Number of columns and lines
+int	nl1, nl2		#I Lines to compute
+int	navg			#I Averaging box half-size
+int	nrej			#I Number of high pixels to reject from average
+int	nbkg			#I Median background width
+real	var0			#I Variance coefficient for DN^0 term
+real	var1			#I Variance coefficient for DN^1 term
+real	var2			#I Variance coefficient for DN^2 term
+int	nsig			#I Sigma box size
+real	lcrsig, hcrsig		#I Threshold sigmas outside of object
+real	lobjsig, hobjsig	#I Object threshold sigmas
+int	crval			#I CR mask value
+int	objval			#I Object mask value
+
+int	i, j, c, c1, c2, c3, c4, l, l1, l2, l3, l4, n1, n2
+int	navg2, nbkg2, nsig2, plsig, phsig
+real	data, avg, bkg, sigma, losig, hosig
+real	low, high, cravg(), amedr()
+pointer	stack, avgs, bkgs, sigs, work1, work2
+pointer	ptr1, ptr2, ip, mp, op, pp, ap, sp
+
+begin
+	navg2 = (2 * navg + 1) ** 2
+	nbkg2 = (2 * (navg + nbkg) + 1) ** 2 - navg2
+	nsig2 = nsig * nsig
+
+	call smark (stack)
+	call salloc (avgs, nc, TY_REAL)
+	call salloc (bkgs, nc, TY_REAL)
+	call salloc (sigs, nc, TY_REAL)
+	call salloc (work1, navg2, TY_REAL)
+	call salloc (work2, max (nsig2, nbkg2), TY_REAL)
+
+	if (var0 != 0. && var1 == 0. && var2 ==0.)
+	    call amovkr (sqrt(var0), Memr[sigs], nc)
+
+	avgs = avgs - 1
+	sigs = sigs - 1
+	bkgs = bkgs - 1
+
+	losig = lobjsig / sqrt (real(navg2-1))
+	hosig = hobjsig / sqrt (real(navg2-1))
+
+	do l = nl1, nl2 {
+	    # Compute statistics.
+	    l1 = max (1, l-navg-nbkg)
+	    l2 = max (1, l-navg)
+	    l3 = min (nl, l+navg)
+	    l4 = min (nl, l+navg+nbkg)
+	    ap = aout + (l - 1) * nc
+	    do c = 1, nc {
+		c1 = max (1, c-navg-nbkg)
+		c2 = max (1, c-navg)
+		c3 = min (nc, c+navg)
+		c4 = min (nc, c+navg+nbkg)
+		ptr1 = work1
+		ptr2 = work2
+		n1 = 0
+		n2 = 0
+		do j = l1, l2-1 {
+		    ip = in + (j - 1) * nc + c1 - 1
+		    mp = pin + (j - 1) * nc + c1 - 1
+		    do i = c1, c4 {
+			if (Mems[mp] == 0) {
+			    Memr[ptr2] = Memr[ip]
+			    n2 = n2 + 1
+			    ptr2 = ptr2 + 1
+			}
+			ip = ip + 1
+			mp = mp + 1
+		    }
+		}
+		do j = l2, l3 {
+		    ip = in + (j - 1) * nc + c1 - 1
+		    mp = pin + (j - 1) * nc + c1 - 1
+		    do i = c1, c2-1 {
+			if (Mems[mp] == 0) {
+			    Memr[ptr2] = Memr[ip]
+			    n2 = n2 + 1
+			    ptr2 = ptr2 + 1
+			}
+			ip = ip + 1
+			mp = mp + 1
+		    }
+		    do i = c2, c3 {
+			if ((j != l || i != c) && Mems[mp] == 0) {
+			    Memr[ptr1] = Memr[ip]
+			    n1 = n1 + 1
+			    ptr1 = ptr1 + 1
+			}
+			ip = ip + 1
+			mp = mp + 1
+		    }
+		    do i = c3+1, c4 {
+			if (Mems[mp] == 0) {
+			    Memr[ptr2] = Memr[ip]
+			    n2 = n2 + 1
+			    ptr2 = ptr2 + 1
+			}
+			ip = ip + 1
+			mp = mp + 1
+		    }
+		}
+		do j = l3+1, l4 {
+		    ip = in + (j - 1) * nc + c1 - 1
+		    mp = pin + (j - 1) * nc + c1 - 1
+		    do i = c1, c4 {
+			if (Mems[mp] == 0) {
+			    Memr[ptr2] = Memr[ip]
+			    n2 = n2 + 1
+			    ptr2 = ptr2 + 1
+			}
+			ip = ip + 1
+		    }
+		}
+		if (n1 > 0)
+		    avg = cravg (Memr[work1], n1, nrej)
+		else
+		    avg = INDEFR
+		if (n2 > 0)
+		    bkg = amedr (Memr[work2], n2)
+		else
+		    bkg = INDEFR
+		Memr[bkgs+c] = bkg
+		Memr[avgs+c] = avg
+		if (aout != NULL) {
+		    if (IS_INDEFR(avg) || IS_INDEFR(bkg))
+			Memr[ap] = 0.
+		    else
+			Memr[ap] = avg - bkg
+		    ap = ap + 1
+		}
+	    }
+
+	    # Compute sigmas and output if desired.
+	    if (var0 != 0. || var1 != 0. || var2 != 0.) {
+		if (var1 != 0.) {
+		    if (var2 != 0.) {
+			do c = 1, nc {
+			    data = max (0., Memr[avgs+c])
+			    Memr[sigs+c] = sqrt (var0+var1*data+var2*data**2)
+			}
+		    } else {
+			do c = 1, nc {
+			    data = max (0., Memr[avgs+c])
+			    Memr[sigs+c] = sqrt (var0 + var1 * data)
+			}
+		    }
+		} else if (var2 != 0.) {
+		    do c = 1, nc {
+			data = max (0., Memr[avgs+c])
+			Memr[sigs+c] = sqrt (var0 + var2 * data**2)
+		    }
+		}
+	    } else {
+		# Compute sigmas from percentiles.  This is done in blocks.
+		if (mod (l-nl1, nsig) == 0 && l<nl-nsig+1) {
+		    do c = 1, nc-nsig+1, nsig {
+			ptr2 = work2
+			n2 = 0
+			do j = l, l+nsig-1 {
+			    ip = in + (j - 1) * nc + c - 1
+			    mp = pin + (j - 1) * nc + c - 1
+			    do i = 1, nsig {
+				if (Mems[mp] == 0) {
+				    Memr[ptr2] = Memr[ip]
+				    n2 = n2 + 1
+				    ptr2 = ptr2 + 1
+				}
+				ip = ip + 1
+				mp = mp + 1
+			    }
+			}
+			if (n2 > 10) {
+			    call asrtr (Memr[work2], Memr[work2], n2)
+			    plsig = nint (PLSIG*n2/100.-1)
+			    phsig = nint (PHSIG*n2/100.-1)
+			    sigma = (Memr[work2+phsig]-Memr[work2+plsig])/2.
+			} else
+			    sigma = INDEFR
+			call amovkr (sigma, Memr[sigs+c], nsig)
+		    }
+		    call amovkr (sigma, Memr[sigs+c], nc-c+1)
+		}
+	    }
+	    if (sout != NULL) {
+		sp = sout + (l - 1) * nc
+		do c = 1, nc {
+		    sigma = Memr[sigs+c]
+		    if (IS_INDEFR(sigma))
+			Memr[sp] = 0.
+		    else
+			Memr[sp] = sigma
+		    sp = sp + 1
+		}
+	    }
+
+	    # Detect, fix, and flag cosmic rays.
+	    if (pout == NULL && out == NULL)
+		;
+	    if (pout == NULL) {
+		ip = in + (l - 1) * nc
+		mp = pin + (l - 1) * nc
+		op = out + (l - 1) * nc
+		do c = 1, nc {
+		    data = Memr[ip]
+		    avg = Memr[avgs+c]
+		    bkg = Memr[bkgs+c]
+		    sigma = Memr[sigs+c]
+		    if (!(Mems[mp] != 0 || IS_INDEFR(avg) ||
+			IS_INDEFR(bkg) || IS_INDEFR(sigma))) {
+			low = bkg - losig * sigma
+			high = bkg + hosig * sigma
+			if (avg < low || avg > high) {
+			    Memr[op] = data
+			} else {
+			    low = avg - lcrsig * sigma
+			    high = avg + hcrsig * sigma
+			    if (data < low || data > high)
+				Memr[op] = avg
+			    else
+				Memr[op] = data
+			}
+		    } else
+			Memr[op] = data
+		    ip = ip + 1
+		    mp = mp + 1
+		    op = op + 1
+		}
+	    } else if (out == NULL) {
+		ip = in + (l - 1) * nc
+		mp = pin + (l - 1) * nc
+		pp = pout + (l - 1) * nc
+		do c = 1, nc {
+		    data = Memr[ip]
+		    avg = Memr[avgs+c]
+		    bkg = Memr[bkgs+c]
+		    sigma = Memr[sigs+c]
+		    if (!(Mems[mp] != 0 || IS_INDEFR(avg) ||
+			IS_INDEFR(bkg) || IS_INDEFR(sigma))) {
+			low = bkg - losig * sigma
+			high = bkg + hosig * sigma
+			if (avg < low || avg > high)
+			    Mems[pp] = objval
+			else {
+			    low = avg - lcrsig * sigma
+			    high = avg + hcrsig * sigma
+			    if (data < low || data > high)
+				Mems[pp] = crval
+			}
+		    }
+		    ip = ip + 1
+		    mp = mp + 1
+		    pp = pp + 1
+		}
+	    } else {
+		ip = in + (l - 1) * nc
+		mp = pin + (l - 1) * nc
+		op = out + (l - 1) * nc
+		pp = pout + (l - 1) * nc
+		do c = 1, nc {
+		    data = Memr[ip]
+		    avg = Memr[avgs+c]
+		    bkg = Memr[bkgs+c]
+		    sigma = Memr[sigs+c]
+		    if (!(Mems[mp] != 0 || IS_INDEFR(avg) ||
+			IS_INDEFR(bkg) || IS_INDEFR(sigma))) {
+			low = bkg - losig * sigma
+			high = bkg + hosig * sigma
+			if (avg < low || avg > high) {
+			    Memr[op] = data
+			    Mems[pp] = objval
+			} else {
+			    low = avg - lcrsig * sigma
+			    high = avg + hcrsig * sigma
+			    if (data < low || data > high) {
+				Memr[op] = avg
+				Mems[pp] = crval
+			    } else
+				Memr[op] = data
+			}
+		    } else
+			Memr[op] = data
+		    ip = ip + 1
+		    mp = mp + 1
+		    op = op + 1
+		    pp = pp + 1
+		}
+	    }
+	}
+
+	call sfree (stack)
+end
+
+
+# CRAVG -- Compute average with the highest nrej points excluded.
+# When nrej is greater than 2 the data array will be returned sorted.
+
+real procedure cravg (data, npts, nrej)
+
+real	data[npts]		#I Input data (will be sorted if nrej>2)
+int	npts			#I Number of data points
+int	nrej			#I Number of data points to reject
+
+int	i
+real	sum, max1, max2, val
+
+begin
+	if (npts <= nrej)
+	    return (INDEFR)
+
+	switch (nrej) {
+	case 0:
+	    sum = 0.
+	    do i = 1, npts
+		sum = sum + data[i]
+	case 1:
+	    sum = 0.
+	    max1 = data[1]
+	    do i = 2, npts {
+		val = data[i]
+		if (val > max1) {
+		    sum = sum + max1
+		    max1 = val
+		} else
+		    sum = sum + val
+	    }
+	case 2:
+	    sum = 0.
+	    max1 = min (data[1], data[2])
+	    max2 = max (data[1], data[2])
+	    do i = 3, npts {
+		val = data[i]
+		if (val > max1) {
+		    sum = sum + max1
+		    if (val > max2) {
+			max1 = max2
+			max2 = val
+		    } else
+			max1 = val
+		} else
+		    sum = sum + val
+	    }
+	default:
+	    call asrtr (data, data, npts)
+	    sum = 0.
+	    do i = 1, npts-nrej
+		sum = sum + data[i]
+	}
+
+	return (sum / (npts - nrej))
+end
diff --git a/noao/imred/crutil/src/t_crgrow.x b/noao/imred/crutil/src/t_crgrow.x
new file mode 100644
index 00000000..7316cc76
--- /dev/null
+++ b/noao/imred/crutil/src/t_crgrow.x
@@ -0,0 +1,182 @@
+include	<error.h>
+include	<imhdr.h>
+
+# T_CRGROW -- Grow cosmic ray mask identifications.
+
+procedure t_crgrow ()
+
+int	input				# Input masks
+int	output				# Output masks
+real	radius				# Radius
+int	inval				# Input mask value to grow
+int	outval				# Output grown mask value
+
+pointer	sp, inmask, outmask, temp1, temp2
+pointer	im, ptr
+
+int	imtopenp(), imtlen(), imtgetim()
+bool	strne()
+int	clgeti()
+real	clgetr()
+pointer	immap()
+errchk	immap, crgrow
+
+begin
+	call smark (sp)
+	call salloc (inmask, SZ_FNAME, TY_CHAR)
+	call salloc (outmask, SZ_FNAME, TY_CHAR)
+	call salloc (temp1, SZ_FNAME, TY_CHAR)
+	call salloc (temp2, SZ_FNAME, TY_CHAR)
+
+	# Task parameters.
+	input = imtopenp ("input")
+	output = imtopenp ("output")
+	radius = max (0., clgetr ("radius"))
+	inval = clgeti ("inval")
+	outval = clgeti ("outval")
+
+	if (imtlen (output) != imtlen (input))
+	    call error (1, "Input and output lists do not match")
+
+	# Grow the cosmic ray masks.
+	while (imtgetim (input, Memc[inmask], SZ_FNAME) != EOF) {
+	    call strcpy (Memc[inmask], Memc[outmask], SZ_FNAME)
+	    if (imtgetim (output, Memc[outmask], SZ_FNAME) == EOF)
+		call error (1, "Output list ended prematurely")
+	    if (strne (Memc[inmask], Memc[outmask])) {
+		call imgcluster (Memc[inmask], Memc[temp1], SZ_FNAME)
+		call imgcluster (Memc[outmask], Memc[temp2], SZ_FNAME)
+		iferr (call imcopy (Memc[temp1], Memc[temp2])) {
+		    call erract (EA_WARN)
+		    next
+		}
+		im = immap (Memc[inmask], READ_ONLY, TY_CHAR)
+		iferr (call imgstr (im, "extname", Memc[temp1], SZ_FNAME))
+		    call strcpy ("pl", Memc[temp1], SZ_FNAME)
+		call imunmap (im)
+	    }
+	    call xt_maskname (Memc[outmask], Memc[temp1], 0, Memc[outmask],
+		SZ_FNAME)
+
+	    if (radius < 1.)
+		next
+
+	    iferr {
+		im = NULL
+	        ptr = immap (Memc[outmask], READ_WRITE, 0); im = ptr
+		call crgrow (im, radius, inval, outval)
+	    } then {
+		call erract (EA_WARN)
+		if (strne (Memc[inmask], Memc[outmask])) {
+		    if (im != NULL) {
+			call imunmap (im)
+			iferr (call imdelete (Memc[outmask]))
+			    call erract (EA_WARN)
+		    }
+		}
+	    }
+
+	    if (im != NULL)
+		call imunmap (im)
+	}
+
+	call imtclose (output)
+	call imtclose (input)
+	call sfree (sp)
+end
+
+
+# CRGROW -- Grow cosmic rays.
+
+procedure crgrow (im, grow, inval, outval)
+
+pointer	im			# Mask pointer (Read/Write)
+real	grow			# Radius (pixels)
+int	inval			# Input mask value for pixels to grow
+int	outval			# Output mask value for grown pixels
+
+int	i, j, k, l, nc, nl, ngrow, nbufs, val1, val2
+long	v1[2], v2[2]
+real	grow2, y2
+pointer	, buf, buf1, buf2, ptr
+
+int	imgnli(), impnli()
+errchk	calloc, imgnli, impnli
+
+begin
+	if (grow < 1. || inval == 0)
+	    return
+
+	grow2 = grow * grow
+	ngrow = int (grow)
+	buf = NULL
+
+	iferr {
+	    if (IM_NDIM(im) > 2)
+		call error (1,
+		    "Only one or two dimensional masks are allowed")
+
+	    nc = IM_LEN(im, 1)
+	    if (IM_NDIM(im) > 1)
+		nl = IM_LEN(im,2)
+	    else
+		nl = 1
+
+	    # Initialize buffering.
+	    nbufs = min (1 + 2 * ngrow, nl)
+	    call calloc (buf, nc*nbufs, TY_INT)
+
+	    call amovkl (long(1), v1, IM_NDIM(im))
+	    call amovkl (long(1), v2, IM_NDIM(im))
+	    while (imgnli (im, buf1, v1) != EOF) {
+		j = v1[2] - 1
+		buf2 = buf + mod (j, nbufs) * nc
+		do i = 1, nc {
+		    val1 = Memi[buf1]
+		    val2 = Memi[buf2]
+		    if ((IS_INDEFI(inval) && val1 != 0) || val1 == inval) {
+			do k = max(1,j-ngrow), min (nl,j+ngrow) {
+			    ptr = buf + mod (k, nbufs) * nc - 1
+			    y2 = (k - j) ** 2
+			    do l = max(1,i-ngrow), min (nc,i+ngrow) {
+				if ((l-i)**2 + y2 > grow2)
+				    next
+				Memi[ptr+l] = -val1
+			    }
+			}
+		    } else {
+			if (val2 >= 0)
+			    Memi[buf2] = val1
+		    }
+		    buf1 = buf1 + 1
+		    buf2 = buf2 + 1
+		}
+
+		if (j > ngrow) {
+		    while (impnli (im, buf2, v2) != EOF) {
+			k = v2[2] - 1
+			buf1 = buf + mod (k, nbufs) * nc
+			do i = 1, nc {
+			    val1 = Memi[buf1]
+			    if (val1 < 0) {
+				if (IS_INDEFI(outval))
+				    Memi[buf2] = -val1
+				else 
+				    Memi[buf2] = outval
+			    } else
+				Memi[buf2] = val1
+			    Memi[buf1] = 0.
+			    buf1 = buf1 + 1
+			    buf2 = buf2 + 1
+			}
+			if (j != nl)
+			    break
+		    }
+		}
+	    }
+	} then
+	    call erract (EA_ERROR)
+
+	if (buf != NULL)
+	    call mfree (buf, TY_INT)
+end
diff --git a/noao/imred/crutil/src/t_crmedian.x b/noao/imred/crutil/src/t_crmedian.x
new file mode 100644
index 00000000..6c7d0fba
--- /dev/null
+++ b/noao/imred/crutil/src/t_crmedian.x
@@ -0,0 +1,417 @@
+include	<imhdr.h>
+include	<mach.h>
+
+define	MAXBUF		500000	# Maximum pixel buffer
+
+define	PLSIG		15.87	# Low percentile
+define	PHSIG		84.13	# High percentile
+
+
+# T_CRMEDIAN -- Detect, fix, and flag cosmic rays.
+# Deviant pixels relative to a local median and sigma are detected and
+# replaced by the median value and/or written to a cosmic ray mask.
+
+procedure t_crmedian ()
+
+pointer	input			# Input image
+pointer	output			# Output image
+pointer	crmask			# Output mask
+pointer	median			# Output median
+pointer	sigma			# Output sigma
+pointer	residual		# Output residual
+real	var0			# Variance coefficient for DN^0 term
+real	var1			# Variance coefficient for DN^1 term
+real	var2			# Variance coefficient for DN^2 term
+real	lsig, hsig		# Threshold sigmas
+int	ncmed, nlmed		# Median box size
+int	ncsig, nlsig		# Sigma box size
+
+int	i, nc, nl, nlstep, l1, l2, l3, l4, nl1
+pointer	sp, extname, in, out, pm, mim, sim, rim
+pointer	inbuf, outbuf, pmbuf, mbuf, sbuf, rbuf
+real	clgetr()
+int	clgeti(), nowhite()
+pointer	immap(), imgs2r(), imps2r(), imps2s()
+errchk	immap, imgs2r, imps2r, imps2s, crmedian, imgstr
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (crmask, SZ_FNAME, TY_CHAR)
+	call salloc (residual, SZ_FNAME, TY_CHAR)
+	call salloc (median, SZ_FNAME, TY_CHAR)
+	call salloc (sigma, SZ_FNAME, TY_CHAR)
+	call salloc (extname, SZ_FNAME, TY_CHAR)
+
+	# Get parameters.
+	call clgstr ("input", Memc[input], SZ_FNAME)
+	call clgstr ("output", Memc[output], SZ_FNAME)
+	call clgstr ("crmask", Memc[crmask], SZ_FNAME)
+	call clgstr ("median", Memc[median], SZ_FNAME)
+	call clgstr ("sigma", Memc[sigma], SZ_FNAME)
+	call clgstr ("residual", Memc[residual], SZ_FNAME)
+	var0 = clgetr ("var0")
+	var1 = clgetr ("var1")
+	var2 = clgetr ("var2")
+	lsig = clgetr ("lsigma")
+	hsig = clgetr ("hsigma")
+	ncmed = clgeti ("ncmed")
+	nlmed = clgeti ("nlmed")
+	ncsig = clgeti ("ncsig")
+	nlsig = clgeti ("nlsig")
+
+	# Map the input and output images.
+	in = NULL; out = NULL; pm = NULL; mim = NULL; sim = NULL; rim = NULL
+	inbuf = NULL; outbuf = NULL; pmbuf = NULL
+	mbuf = NULL; sbuf=NULL; rbuf = NULL
+	in = immap (Memc[input], READ_ONLY, 0)
+	if (nowhite (Memc[output], Memc[output], SZ_FNAME) > 0)
+	    out = immap (Memc[output], NEW_COPY, in)
+	if (nowhite (Memc[crmask], Memc[crmask], SZ_FNAME) > 0) {
+	    if (Memc[crmask] == '!')
+		call imgstr (in, Memc[crmask+1], Memc[crmask], SZ_FNAME)
+	    iferr (call imgstr (in, "extname", Memc[extname], SZ_FNAME))
+		call strcpy ("pl", Memc[extname], SZ_FNAME)
+	    call xt_maskname (Memc[crmask], Memc[extname], 0, Memc[crmask],
+		SZ_FNAME)
+	    pm = immap (Memc[crmask], NEW_COPY, in)
+	}
+	if (nowhite (Memc[median], Memc[median], SZ_FNAME) > 0)
+	    mim = immap (Memc[median], NEW_COPY, in)
+	if (nowhite (Memc[sigma], Memc[sigma], SZ_FNAME) > 0)
+	    sim = immap (Memc[sigma], NEW_COPY, in)
+	if (nowhite (Memc[residual], Memc[residual], SZ_FNAME) > 0)
+	    rim = immap (Memc[residual], NEW_COPY, in)
+
+	# Go through the input in large blocks of lines.  If the
+	# block is smaller than the whole image overlap the blocks
+	# so the median only has boundaries at the ends of the image.
+	# However, the output is done in non-overlapping blocks with
+	# the pointers are adjusted so that addresses can be in the space
+	# of the input block.  CRMEDIAN does not address outside of
+	# the output data block.  Set the mask values based on the
+	# distances to the nearest good pixels.
+
+	nc = IM_LEN(in,1)
+	nl = IM_LEN(in,2)
+	nlstep = max (1, MAXBUF / nc - nlmed)
+
+	do i = 1, nl, nlstep {
+	    l1 = i
+	    l2 = min (nl, i + nlstep - 1)
+	    l3 = max (1, l1 - nlmed / 2)
+	    l4 = min (nl, l2 + nlmed / 2)
+	    nl1 = l4 - l3 + 1
+	    inbuf = imgs2r (in, 1, nc, l3, l4)
+	    if (out != NULL)
+		outbuf = imps2r (out, 1, nc, l1, l2) - (l1 - l3) * nc
+	    if (pm != NULL)
+		pmbuf = imps2s (pm, 1, nc, l1, l2) - (l1 - l3) * nc
+	    if (mim != NULL)
+		mbuf = imps2r (mim, 1, nc, l1, l2) - (l1 - l3) * nc
+	    if (sim != NULL)
+		sbuf = imps2r (sim, 1, nc, l1, l2) - (l1 - l3) * nc
+	    if (rim != NULL)
+		rbuf = imps2r (rim, 1, nc, l1, l2) - (l1 - l3) * nc
+	    call crmedian (inbuf, outbuf, pmbuf, mbuf, sbuf, rbuf,
+		nc, nl1, l1-l3+1, l2-l3+1, ncmed, nlmed, var0, var1, var2,
+		ncsig, nlsig, lsig, hsig)
+	}
+
+	if (rim != NULL)
+	    call imunmap (rim)
+	if (sim != NULL)
+	    call imunmap (sim)
+	if (mim != NULL)
+	    call imunmap (mim)
+	if (pm != NULL)
+	    call imunmap (pm)
+	if (out != NULL)
+	    call imunmap (out)
+	call imunmap (in)
+	call sfree (sp)
+end
+
+
+# CRMEDIAN -- Detect, replace, and flag cosmic rays.
+# A local background is computed using moving box medians to avoid
+# contaminating bad pixels.  If variance model is given then that is
+# used otherwise a local sigma is computed in blocks (it is not a moving box
+# for efficiency) by using a percentile point of the sorted pixel values to
+# estimate the width of the distribution uncontaminated by bad pixels).  Once
+# the background and sigma are known deviant pixels are found by using sigma
+# threshold factors.
+
+procedure crmedian (in, out, pout, mout, sout, rout, nc, nl, nl1, nl2,
+	ncmed, nlmed, var0, var1, var2, ncsig, nlsig, lsig, hsig)
+
+pointer	in			#I Input data
+pointer	out			#O Output data
+pointer	pout			#O Output mask (0=good, 1=bad)
+pointer	mout			#O Output medians
+pointer	sout			#O Output sigmas
+pointer	rout			#O Output residuals
+int	nc, nl			#I Number of columns and lines
+int	nl1, nl2		#I Lines to compute
+int	ncmed, nlmed		#I Median box size
+real	var0			# Variance coefficient for DN^0 term
+real	var1			# Variance coefficient for DN^1 term
+real	var2			# Variance coefficient for DN^2 term
+int	ncsig, nlsig		#I Sigma box size
+real	lsig, hsig		#I Threshold sigmas
+
+int	i, j, k, l, m, plsig, phsig
+real	data, med, sigma, low, high, amedr()
+pointer	stack, meds, sigs, work, ptr, ip, op, pp, mp, sp, rp
+
+begin
+	call smark (stack)
+	call salloc (meds, nc, TY_REAL)
+	call salloc (sigs, nc, TY_REAL)
+	call salloc (work, max (ncsig*nlsig, ncmed*nlmed), TY_REAL)
+
+	if (var0 != 0. && var1 == 0. && var2 ==0.)
+	    call amovkr (sqrt(var0), Memr[sigs], nc)
+
+	meds = meds - 1
+	sigs = sigs - 1
+
+	i = ncsig * nlsig
+	plsig = nint (PLSIG*i/100.-1)
+	phsig = nint (PHSIG*i/100.-1)
+
+	do i = nl1, nl2 {
+
+	    # Compute median and output if desired.  This is a moving median.
+	    l = min (nl, i+nlmed/2)
+	    l = max (1, l-nlmed+1)
+	    mp = mout + (i - 1) * nc 
+	    do j = 1, nc {
+		k = min (nc, j+ncmed/2)
+		k = max (1, k-ncmed+1)
+		ptr = work
+		ip = in + (l - 1) * nc + k - 1
+		do m = l, l+nlmed-1 {
+		    call amovr (Memr[ip], Memr[ptr], ncmed)
+		    ip = ip + nc
+		    ptr = ptr + ncmed
+		}
+		med = amedr (Memr[work], ncmed * nlmed)
+		Memr[meds+j] = med
+		if (mout != NULL) {
+		    Memr[mp] = med
+		    mp = mp + 1
+		}
+	    }
+
+	    # Compute sigmas and output if desired.
+	    if (var0 != 0. || var1 != 0. || var2 != 0.) {
+		if (var1 != 0.) {
+		    if (var2 != 0.) {
+			do j = 1, nc {
+			    data = max (0., Memr[meds+j])
+			    Memr[sigs+j] = sqrt (var0 + var1*data + var2*data**2)
+			}
+		    } else {
+			do j = 1, nc {
+			    data = max (0., Memr[meds+j])
+			    Memr[sigs+j] = sqrt (var0 + var1 * data)
+			}
+		    }
+		} else if (var2 != 0.) {
+		    do j = 1, nc {
+			data = max (0., Memr[meds+j])
+			Memr[sigs+j] = sqrt (var0 + var2 * data**2)
+		    }
+		}
+	    } else {
+		# Compute sigmas from percentiles.  This is done in blocks.
+		if (mod (i-nl1, nlsig) == 0 && i<nl-nlsig+1) {
+		    do j = 1, nc-ncsig+1, ncsig {
+			ptr = work
+			ip = in + (i - 1) * nc + j - 1
+			do k = i, i+nlsig-1 {
+			    call amovr (Memr[ip], Memr[ptr], ncsig)
+			    ip = ip + nc
+			    ptr = ptr + ncsig
+			}
+			call asrtr (Memr[work], Memr[work], ncsig*nlsig)
+			sigma = (Memr[work+phsig] - Memr[work+plsig]) / 2.
+			call amovkr (sigma, Memr[sigs+j], ncsig)
+		    }
+		    call amovkr (sigma, Memr[sigs+j], nc-j+1)
+		}
+	    }
+	    if (sout != NULL) {
+		sp = sout + (i - 1) * nc
+		do j = 1, nc {
+		    Memr[sp] = Memr[sigs+j]
+		    sp = sp + 1
+		}
+	    }
+
+	    # Detect, fix, and flag cosmic rays.
+	    if (rout == NULL) {
+		if (pout == NULL) {
+		    ip = in + (i - 1) * nc
+		    op = out + (i - 1) * nc
+		    do j = 1, nc {
+			data = Memr[ip]
+			med = Memr[meds+j]
+			sigma = Memr[sigs+j]
+			low = med - lsig * sigma
+			high = med + hsig * sigma
+			if (data < low || data > high)
+			    Memr[op] = med
+			else
+			    Memr[op] = data
+			ip = ip + 1
+			op = op + 1
+		    }
+		} else if (out == NULL) {
+		    ip = in + (i - 1) * nc
+		    pp = pout + (i - 1) * nc
+		    do j = 1, nc {
+			data = Memr[ip]
+			med = Memr[meds+j]
+			sigma = Memr[sigs+j]
+			low = med - lsig * sigma
+			high = med + hsig * sigma
+			if (data < low || data > high)
+			    Mems[pp] = 1
+			else
+			    Mems[pp] = 0
+			ip = ip + 1
+			pp = pp + 1
+		    }
+		} else {
+		    ip = in + (i - 1) * nc
+		    op = out + (i - 1) * nc
+		    pp = pout + (i - 1) * nc
+		    do j = 1, nc {
+			data = Memr[ip]
+			med = Memr[meds+j]
+			sigma = Memr[sigs+j]
+			low = med - lsig * sigma
+			high = med + hsig * sigma
+			if (data < low || data > high) {
+			    Memr[op] = med
+			    Mems[pp] = 1
+			} else {
+			    Memr[op] = data
+			    Mems[pp] = 0
+			}
+			ip = ip + 1
+			op = op + 1
+			pp = pp + 1
+		    }
+		}
+	    } else {
+		if (pout == NULL && out == NULL) {
+		    ip = in + (i - 1) * nc
+		    rp = rout + (i - 1) * nc
+		    do j = 1, nc {
+			data = Memr[ip]
+			med = Memr[meds+j]
+			sigma = Memr[sigs+j]
+			if (sigma > 0.)
+			    Memr[rp] = (data - med) / sigma
+			else {
+			    if ((data - med) < 0.)
+				Memr[rp] = -MAX_REAL
+			    else
+				Memr[rp] = MAX_REAL
+			}
+			ip = ip + 1
+			rp = rp + 1
+		    }
+		} else if (pout == NULL) {
+		    ip = in + (i - 1) * nc
+		    op = out + (i - 1) * nc
+		    rp = rout + (i - 1) * nc
+		    do j = 1, nc {
+			data = Memr[ip]
+			med = Memr[meds+j]
+			sigma = Memr[sigs+j]
+			low = med - lsig * sigma
+			high = med + hsig * sigma
+			if (data < low || data > high)
+			    Memr[op] = med
+			else
+			    Memr[op] = data
+			if (sigma > 0.)
+			    Memr[rp] = (data - med) / sigma
+			else {
+			    if ((data - med) < 0.)
+				Memr[rp] = -MAX_REAL
+			    else
+				Memr[rp] = MAX_REAL
+			}
+			ip = ip + 1
+			op = op + 1
+			rp = rp + 1
+		    }
+		} else if (out == NULL) {
+		    ip = in + (i - 1) * nc
+		    pp = pout + (i - 1) * nc
+		    rp = rout + (i - 1) * nc
+		    do j = 1, nc {
+			data = Memr[ip]
+			med = Memr[meds+j]
+			sigma = Memr[sigs+j]
+			low = med - lsig * sigma
+			high = med + hsig * sigma
+			if (data < low || data > high)
+			    Mems[pp] = 1
+			else
+			    Mems[pp] = 0
+			if (sigma > 0.)
+			    Memr[rp] = (data - med) / sigma
+			else {
+			    if ((data - med) < 0.)
+				Memr[rp] = -MAX_REAL
+			    else
+				Memr[rp] = MAX_REAL
+			}
+			ip = ip + 1
+			pp = pp + 1
+			rp = rp + 1
+		    }
+		} else {
+		    ip = in + (i - 1) * nc
+		    op = out + (i - 1) * nc
+		    pp = pout + (i - 1) * nc
+		    rp = rout + (i - 1) * nc
+		    do j = 1, nc {
+			data = Memr[ip]
+			med = Memr[meds+j]
+			sigma = Memr[sigs+j]
+			low = med - lsig * sigma
+			high = med + hsig * sigma
+			if (data < low || data > high) {
+			    Memr[op] = med
+			    Mems[pp] = 1
+			} else {
+			    Memr[op] = data
+			    Mems[pp] = 0
+			}
+			if (sigma > 0.)
+			    Memr[rp] = (data - med) / sigma
+			else {
+			    if ((data - med) < 0.)
+				Memr[rp] = -MAX_REAL
+			    else
+				Memr[rp] = MAX_REAL
+			}
+			ip = ip + 1
+			op = op + 1
+			pp = pp + 1
+			rp = rp + 1
+		    }
+		}
+	    }
+	}
+
+	call sfree (stack)
+end
diff --git a/noao/imred/crutil/src/x_crutil.x b/noao/imred/crutil/src/x_crutil.x
new file mode 100644
index 00000000..48408a9f
--- /dev/null
+++ b/noao/imred/crutil/src/x_crutil.x
@@ -0,0 +1,4 @@
+task	cosmicrays	= t_cosmicrays,
+	craverage	= t_craverage,
+	crgrow		= t_crgrow,
+	crmedian	= t_crmedian
diff --git a/noao/imred/crutil/src/xtmaskname.x b/noao/imred/crutil/src/xtmaskname.x
new file mode 100644
index 00000000..eb132e12
--- /dev/null
+++ b/noao/imred/crutil/src/xtmaskname.x
@@ -0,0 +1,125 @@
+# MASKNAME -- Make a mask name.  This creates a FITS mask extension if
+# possible, otherwise it creates a pixel list file.  To create a FITS
+# extension the filename must explicitly select the FITS kernel or the
+# default image type must be a FITS file.  The input and output strings
+# may be the same.
+
+procedure xt_maskname (fname, extname, mode, mname, maxchar)
+
+char	fname[ARB]			#I File name
+char	extname[ARB]			#I Default pixel mask extension name
+int	mode				#I Mode
+char	mname[maxchar]			#O Output mask name
+int	maxchar				#I Maximum characters in mask name
+
+int	i, fits
+pointer	sp, temp
+
+bool	streq()
+int	strmatch(), stridxs(), strldxs(), strncmp()
+int	envfind(), access(), imaccess()
+
+begin
+	call smark (sp)
+	call salloc (temp, maxchar, TY_CHAR)
+
+	# Determine whether to use FITS pixel mask extensions.  One may set
+	# fits=NO to force use of pl even when FITS mask extensions are
+	# supported.
+	fits = YES
+	if (fits == YES && envfind ("masktype", Memc[temp], maxchar) > 0) {
+	    if (streq (Memc[temp], "pl"))
+		fits = NO
+	}
+	i = strldxs ("]", fname)
+
+	# Check for explicit .pl extension.
+	if (strmatch (fname, ".pl$") > 0)
+	    call strcpy (fname, mname, maxchar)
+
+	# Check for explicit mask extension.
+	else if (strmatch (fname, "type=mask") > 0)
+	    call strcpy (fname, mname, maxchar)
+	else if (strmatch (fname, "type\\\=mask") > 0)
+	    call strcpy (fname, mname, maxchar)
+
+	# Add mask type.
+	else if (i > 0) {
+	    if (mode != READ_ONLY) {
+		call strcpy (fname[i], Memc[temp], maxchar)
+		call sprintf (mname[i], maxchar-i, ",type=mask%s")
+		    call pargstr (Memc[temp])
+	    }
+	    i = stridxs ("[", mname)
+	    if (i > 0) {
+		call strcpy (mname[i+1], Memc[temp], SZ_FNAME)
+		if (strncmp (mname[i+1], "append", 6)==0 ||
+		    strncmp (mname[i+1], "inherit", 7)==0) {
+		    call sprintf (mname[i+1], SZ_FNAME-i+1, "%s,%s")
+			call pargstr (extname)
+			call pargstr (Memc[temp])
+		}
+	    }
+
+	# Create output from rootname name.
+	} else if (fits == YES) {
+	    call strcpy (fname, Memc[temp], SZ_FNAME)
+	    if (mode == READ_ONLY) {
+		call sprintf (mname, maxchar, "%s[%s]")
+		    call pargstr (Memc[temp])
+		    call pargstr (extname)
+	    } else {
+		call sprintf (mname, maxchar, "%s[%s,type=mask]")
+		    call pargstr (Memc[temp])
+		    call pargstr (extname)
+	    }
+	} else
+	    call strcat (".pl", mname, maxchar)
+
+	# Check if extension name is actually given.
+	i = stridxs ("[", mname)
+	if (i > 0) {
+	    call strcpy (mname[i+1], Memc[temp], SZ_FNAME)
+	    if (strncmp (mname[i+1], "append", 6)==0 ||
+		strncmp (mname[i+1], "inherit", 7)==0) {
+		call sprintf (mname[i+1], SZ_FNAME-i+1, "%s,%s")
+		    call pargstr (extname)
+		    call pargstr (Memc[temp])
+	    }
+	}
+		
+	# Convert to pl form if required.
+	i = stridxs ("[", mname)
+	if (i > 0 && mode == READ_ONLY)
+	    fits = imaccess (mname, mode)
+	if (fits == NO && i > 0) {
+	    mname[i] = EOS
+	    if (mode == NEW_IMAGE) {
+		if (access (mname, 0, 0) == NO) {
+		    ifnoerr (call fmkdir (mname))
+			mname[i] =  '/'
+		    else
+			mname[i] = '.'
+		} else
+		    mname[i] =  '/'
+	    } else {
+		if (access (mname, 0, 0) == NO)
+		    mname[i] = '.'
+		else
+		    mname[i] =  '/'
+	    }
+	        
+	    if (strncmp (mname[i+1], "type", 4) == 0 ||
+	        strncmp (mname[i+1], "append", 6) == 0 ||
+	        strncmp (mname[i+1], "inherit", 7) == 0) {
+		mname[i+1] = EOS
+		call strcat (extname, mname, maxchar)
+	    } else {
+		i = stridxs (",]", mname)
+		mname[i] = EOS
+	    }
+	    call strcat (".pl", mname, maxchar)
+	}
+
+	call sfree (sp)
+end