Repatch (from linux) of OSX IRAF

9 files changed, 3872 insertions, 0 deletions

diff --git a/pkg/xtools/fixpix/mkpkg b/pkg/xtools/fixpix/mkpkg
new file mode 100644
index 00000000..4d91ae71
--- /dev/null
+++ b/pkg/xtools/fixpix/mkpkg
@@ -0,0 +1,25 @@
+# XT_FIXPIX package.
+
+$checkout libxtools.a lib$
+$update	  libxtools.a
+$checkin  libxtools.a lib$
+$exit
+
+generic:
+	$set	GEN = "$$generic -k"
+	$ifolder (xtfp.x, xtfp.gx)
+	    $(GEN) xtfp.gx -o xtfp.x $endif
+	;
+
+libxtools.a:
+	$ifeq (USE_GENERIC, yes) $call generic $endif
+
+	setfp.x		<imhdr.h> <imset.h> <pmset.h>
+	xtfixpix.x	<imhdr.h> <imset.h> <pmset.h> xtfixpix.h
+	xtfp.x		<imhdr.h> <pmset.h> xtfixpix.h
+	xtpmmap.x	<ctype.h> <error.h> <imhdr.h> <imset.h> <mach.h>\
+			<mwset.h> <pmset.h>
+	ytfixpix.x	<imhdr.h> <imset.h> <pmset.h> xtfixpix.h
+	ytpmmap.x	<ctype.h> <error.h> <imhdr.h> <imset.h> <mach.h>\
+			<mwset.h> <pmset.h>
+	;
diff --git a/pkg/xtools/fixpix/setfp.x b/pkg/xtools/fixpix/setfp.x
new file mode 100644
index 00000000..5fe2f5c1
--- /dev/null
+++ b/pkg/xtools/fixpix/setfp.x
@@ -0,0 +1,72 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<imset.h>
+include	<pmset.h>
+
+
+# SET_FP -- Set the fixpix mask.
+#
+# This routine transforms the input mask values into the output mask
+# values.  It allows the input mask to have two classes of bad pixels;
+# those which are interpolated and those which are not.
+
+procedure set_fp (im, fp)
+
+pointer	im			#I Input mask image pointer
+pointer	fp			#O FIXPIX interpolation pointer
+
+int	i, j, nc, nl
+long	v[2]
+pointer	data1, data2, pm, pmi
+
+int	imstati(), pm_newcopy()
+pointer	yt_fpinit()
+errchk	malloc, yt_fpinit
+
+begin
+	# Set the image size and data buffers.
+	nc = IM_LEN(im,1)
+	nl = IM_LEN(im,2)
+	call malloc (data1, nc, TY_SHORT)
+	call malloc (data2, nc, TY_SHORT)
+
+	# Get the pixel mask from the image.
+	pm = imstati (im, IM_PMDES)
+
+	# Extract the pixels to be interpolated.
+	pmi = pm_newcopy (pm)
+	v[1] = 1
+	do j = 1, nl {
+	    v[2] = j
+	    call pmglps (pm, v, Mems[data1], 0, nc, PIX_SRC)
+	    do i = 0, nc-1 {
+	        if (Mems[data1+i] > 1)
+		    Mems[data1+i] = 0
+	    }
+	    call pmplps (pmi, v, Mems[data1], 0, nc, PIX_SRC)
+	}
+
+	# Set the interpolation.
+	fp = yt_fpinit (pmi, 2, 3)
+
+	# Merge back the bad pixels which are not interpolated.
+	v[1] = 1
+	do j = 1, nl {
+	    v[2] = j
+	    call pmglps (pm, v, Mems[data1], 0, nc, PIX_SRC)
+	    call pmglps (pmi, v, Mems[data2], 0, nc, PIX_SRC)
+	    do i = 0, nc-1 {
+	        if (Mems[data2+i] != 0)
+		    Mems[data1+i] = Mems[data2+i]
+		else if (Mems[data1+i] > 1)
+		    Mems[data1+i] = 6
+	    }
+	    call pmplps (pm, v, Mems[data1], 0, nc, PIX_SRC)
+	}
+
+	# Finish up.
+	call mfree (data1, TY_SHORT)
+	call mfree (data2, TY_SHORT)
+	#call pm_close (pmi)
+end
diff --git a/pkg/xtools/fixpix/xtfixpix.h b/pkg/xtools/fixpix/xtfixpix.h
new file mode 100644
index 00000000..de30f65d
--- /dev/null
+++ b/pkg/xtools/fixpix/xtfixpix.h
@@ -0,0 +1,24 @@
+# XT_FIXPIX data structure.
+define	FP_LEN		13		# Length of FP structure
+define	FP_PM		Memi[$1]	# Pixel mask pointer
+define	FP_LVAL		Memi[$1+1]	# Mask value for line interpolation
+define	FP_CVAL		Memi[$1+2]	# Mask value for column interpolation
+define	FP_NCOLS	Memi[$1+3]	# Number of columns to interpolate
+define	FP_PCOL		Memi[$1+4]	# Pointer to columns
+define	FP_PL1		Memi[$1+5]	# Pointer to start lines
+define	FP_PL2		Memi[$1+6]	# Pointer to end lines
+define	FP_PV1		Memi[$1+7]	# Pointer to start values	
+define	FP_PV2		Memi[$1+8]	# Pointer to end values	
+define	FP_LMIN		Memi[$1+9]	# Minimum line
+define	FP_LMAX		Memi[$1+10]	# Maximum line
+define	FP_PIXTYPE	Memi[$1+11]	# Pixel type for values
+define	FP_DATA		Memi[$1+12]	# Data values
+
+define	FP_COL		Memi[FP_PCOL($1)+$2-1]
+define	FP_L1		Memi[FP_PL1($1)+$2-1]
+define	FP_L2		Memi[FP_PL2($1)+$2-1]
+define	FP_V1		(FP_PV1($1)+$2-1)
+define	FP_V2		(FP_PV2($1)+$2-1)
+
+define	FP_LDEF		1		# Default line interpolation code
+define	FP_CDEF		2		# Default column interpolation code
diff --git a/pkg/xtools/fixpix/xtfixpix.x b/pkg/xtools/fixpix/xtfixpix.x
new file mode 100644
index 00000000..500824b5
--- /dev/null
+++ b/pkg/xtools/fixpix/xtfixpix.x
@@ -0,0 +1,270 @@
+include	<imhdr.h>
+include	<imset.h>
+include	<pmset.h>
+include	"xtfixpix.h"
+
+
+# XT_FPINIT -- Initialize FIXPIX data structure.
+# If the mask is null or empty a null pointer is returned.
+# If the mask is not empty the mask is examined for bad pixels requiring
+# column interpolation.  The columns and interpolation endpoints are
+# recorded.  Note that line interpolation does not need to be mapped since
+# this can be done efficiently as the reference image is accessed line by
+# line.
+
+pointer procedure xt_fpinit (pm, lvalin, cvalin)
+
+pointer	pm			#I Pixel mask
+int	lvalin			#I Input line interpolation code
+int	cvalin			#I Input column interpolation code
+
+int	i, j, k, l, n, nc, nl, l1, l2, lmin, lmax, ncols, lval, cval, ncompress
+short	val
+long	v[IM_MAXDIM]
+pointer	fp, ptr, col, pl1, pl2
+pointer	sp, buf, cols
+
+bool	pm_empty()
+errchk	pmglrs, pmplrs()
+
+begin
+	# Check for empty mask.
+	if (pm == NULL)
+	    return (NULL)
+	if (pm_empty (pm))
+	    return (NULL)
+
+	# Get mask size.
+	call pm_gsize (pm, i, v, j)
+	nc = v[1]
+	nl = v[2]
+
+	# Allocate memory and data structure.
+	call smark (sp)
+	call salloc (buf, 3*max(nc, nl), TY_SHORT) 
+	call salloc (cols, nc, TY_SHORT)
+	call calloc (fp, FP_LEN, TY_STRUCT)
+
+	# Set the mask codes.  Go through the mask and change any mask codes
+	# that match the input mask code to the output mask code (if they are
+	# different).  This is done to move the mask codes to a range that
+	# won't conflict with the length values.  For any other code replace
+	# the value by the length of the bad region along the line.  This
+	# value will be used in comparison to the length along the column for
+	# setting the interpolation for the narrower dimension.
+
+	if ((IS_INDEFI(lvalin)||lvalin<1) && (IS_INDEFI(cvalin)||cvalin<1)) {
+	    lval = FP_LDEF
+	    cval = FP_CDEF
+	} else if (IS_INDEFI(lvalin) || lvalin < 1) {
+	    lval = FP_LDEF
+	    cval = mod (cvalin - 1, nc) + 1
+	    if (lval == cval)
+		lval = FP_CDEF
+	} else if (IS_INDEFI(cvalin) || cvalin < 1) {
+	    lval = mod (lvalin - 1, nc) + 1
+	    cval = FP_CDEF
+	    if (cval == lval)
+		cval = FP_LDEF
+	} else if (lvalin != cvalin) {
+	    lval = mod (lvalin - 1, nc) + 1
+	    cval = mod (cvalin - 1, nc) + 1
+	} else {
+	    call mfree (fp, TY_STRUCT)
+	    call sfree (sp)
+	    call error (1, "Interpolation codes cannot be the same")
+	}
+	call xt_fpsinterp (pm, nc, nl, v, Mems[buf], lvalin, cvalin, lval, cval)
+
+	# Go through and check if there is any need for column interpolation;
+	# i.e. are there any mask values different from the line interpolation.
+
+	call aclrs (Mems[cols], nc)
+	call amovkl (long(1), v, IM_MAXDIM)
+	do l = 1, nl {
+	    v[2] = l
+	    call pmglrs (pm, v, Mems[buf], 0, nc, 0)
+	    ptr = buf + 3
+	    do i = 2, Mems[buf] {
+		val = Mems[ptr+2]
+		if (val != lval) {
+		    val = 1
+		    n = Mems[ptr+1] 
+		    call amovks (val, Mems[cols+Mems[ptr]-1], n)
+		}
+		ptr = ptr + 3
+	    }
+	}
+	n = 0
+	do i = 1, nc
+	    if (Mems[cols+i-1] != 0)
+		n = n + 1
+
+	# If there are mask codes for either column interpolation or
+	# interpolation lengths along lines to compare against column
+	# interpolation check the interpolation length against the
+	# column and set the line interpolation endpoints to use.
+	# compute the minimum and maximum lines that are endpoints
+	# to restrict the random access pass that will be needed to
+	# get the endpoint values.
+
+	if (n > 0) {
+	    n = n + 10
+	    call malloc (col, n, TY_INT)
+	    call malloc (pl1, n, TY_INT)
+	    call malloc (pl2, n, TY_INT)
+	    ncols = 0
+	    lmin = nl
+	    lmax = 0
+	    ncompress = 0
+	    do i = 1, nc {
+		if (Mems[cols+i-1] == 0)
+		    next
+		v[1] = i
+		do l = 1, nl {
+		    v[2] = l
+		    call pmglps (pm, v, Mems[buf+l-1], 0, 1, 0)
+		}
+		for (l1=1; l1<=nl && Mems[buf+l1-1]==0; l1=l1+1)
+		    ;
+		while (l1 <= nl) {
+		    l1 = l1 - 1
+		    for (l2=l1+1; l2<=nl && Mems[buf+l2-1]!=0; l2=l2+1)
+			;
+		    j = 0
+		    k = nc + l2 - l1 - 1
+		    do l = l1+1, l2-1 {
+			val = Mems[buf+l-1]
+			if (val == cval)
+			    j = j + 1
+			else if (val > nc) {
+			    if (val > k) {
+				j = j + 1
+				val = cval
+			    } else
+				val = lval
+			    v[2] = l
+			    call pmplps (pm, v, val, 0, 1, PIX_SRC)
+			    ncompress = ncompress + 1
+			}
+		    }
+		    if (ncompress > 100) {
+			call pm_compress (pm)
+			ncompress = 0
+		    }
+		    if (j > 0) {
+			if (ncols == n) {
+			    n = n + 10
+			    call realloc (col, n, TY_INT)
+			    call realloc (pl1, n, TY_INT)
+			    call realloc (pl2, n, TY_INT)
+			}
+			j = 1 + l1 - 1
+			k = 1 + l2 - 1
+			lmin = min (lmin, j, k)
+			lmax = max (lmax, j, k)
+			Memi[col+ncols] = i
+			Memi[pl1+ncols] = j
+			Memi[pl2+ncols] = k
+			ncols = ncols + 1
+		    }
+		    for (l1=l2+1; l1<=nl && Mems[buf+l1-1]==0; l1=l1+1)
+			;
+		}
+	    }
+
+	    FP_LMIN(fp) = lmin
+	    FP_LMAX(fp) = lmax
+	    FP_NCOLS(fp) = ncols
+	    FP_PCOL(fp) = col
+	    FP_PL1(fp) = pl1
+	    FP_PL2(fp) = pl2
+	}
+
+	FP_PM(fp) = pm
+	FP_LVAL(fp) = lval
+	FP_CVAL(fp) = cval
+
+	call sfree (sp)
+	return (fp)
+end
+
+
+# XT_SINTERP -- Set length of line interpolation regions.
+# The mask values are set to the length of any column interpolation
+# plus an offset leaving any line and column interpolation codes
+# unchanged.  These values will be used in a second pass to compare
+# to the lengths of line interpolation and then the mask values will
+# be reset to one of the line or column interpolation codes based on
+# the minimum distance.
+
+procedure xt_fpsinterp (pm, nc, nl, v, data, lvalin, cvalin, lvalout, cvalout)
+
+pointer	pm		#I Pixel mask
+int	nc, nl		#I Mask size
+long	v[ARB]		#I Coordinate vector
+short	data[ARB]	#I Data buffer
+int	lvalin		#I Input line interpolation code
+int	cvalin		#I Input column interpolation code
+int	lvalout		#I Output line interpolation code
+int	cvalout		#I Output column interpolation code
+
+int	c, l, c1, c2, val
+bool	pm_linenotempty()
+
+begin
+	call amovkl (long(1), v, IM_MAXDIM)
+	do l = 1, nl {
+	    v[2] = l
+	    if (!pm_linenotempty (pm, v))
+		next
+	    
+	    call pmglps (pm, v, data, 0, nc, 0)
+
+	    for (c1=1; c1<=nc && data[c1]==0; c1=c1+1)
+		;
+	    while (c1 <= nc) {
+		for (c2=c1+1; c2<=nc && data[c2]!=0; c2=c2+1)
+		    ;
+		c2 = c2 - 1
+		do c = c1, c2 {
+		    val = data[c]
+		    if (val == lvalin) {
+			if (lvalin != lvalout)
+			    data[c] = lvalout
+		    } else if (val == cvalin) {
+			if (cvalin != cvalout)
+			    data[c] = cvalout
+		    } else {
+			data[c] = nc + c2 - c1 + 1
+		    }
+		}
+		for (c1=c2+2; c1<=nc && data[c1]==0; c1=c1+1)
+		    ;
+	    }
+
+	    call pmplps (pm, v, data, 0, nc, PIX_SRC)
+	}
+end
+
+
+# XT_FPFREE -- Free FIXPIX data structures.
+
+procedure xt_fpfree (fp)
+
+pointer	fp		#U FIXPIX data structure
+
+begin
+	if (fp == NULL)
+	    return
+	call mfree (FP_PCOL(fp), TY_INT)
+	call mfree (FP_PL1(fp), TY_INT)
+	call mfree (FP_PL2(fp), TY_INT)
+	if (FP_PV1(fp) != NULL)
+	    call mfree (FP_PV1(fp), FP_PIXTYPE(fp))
+	if (FP_PV2(fp) != NULL)
+	    call mfree (FP_PV2(fp), FP_PIXTYPE(fp))
+	if (FP_DATA(fp) != NULL)
+	    call mfree (FP_DATA(fp), FP_PIXTYPE(fp))
+	call mfree (fp, TY_STRUCT)
+end
diff --git a/pkg/xtools/fixpix/xtfp.gx b/pkg/xtools/fixpix/xtfp.gx
new file mode 100644
index 00000000..70893ff8
--- /dev/null
+++ b/pkg/xtools/fixpix/xtfp.gx
@@ -0,0 +1,275 @@
+include	<imhdr.h>
+include	<pmset.h>
+include	"xtfixpix.h"
+
+
+$for (silrd)
+
+# XT_FP -- Get the specified line of image data and replace bad pixels by
+# interpolation.
+
+pointer procedure xt_fp$t (fp, im, line, fd)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+int	fd			#I File descriptor for pixel list
+
+int	col1, col2		#I Section of interest
+int	line1, line2		#I Section of interest
+
+pointer	imgl2$t(), xt_fps$t()
+
+begin
+	# If there are no bad pixels just get the image line and return.
+	if (fp == NULL)
+	    return (imgl2$t (im, line))
+
+	col1 = 1
+	col2 = IM_LEN(im,1)
+	line1 = 1
+	line2 = IM_LEN(im,2)
+
+	return (xt_fps$t (fp, im, line, col1, col2, line1, line2, fd))
+end
+
+
+# XT_FXS -- Get the specified line of image data and replace bad pixels by
+# interpolation within a specified section.
+
+pointer procedure xt_fps$t (fp, im, line, col1, col2, line1, line2, fd)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+int	fd			#I File descriptor for pixel list
+
+int	col1, col2		#I Section of interest
+int	line1, line2		#I Section of interest
+
+int	i, j, nc, nl, ncols, c1, c2, l1, l2, l3, l4
+long	v[IM_MAXDIM]
+$if (datatype == silr)
+real	a, b, c, d, val
+$else
+PIXEL	a, b, c, d, val
+$endif
+PIXEL	indef
+pointer	pm, data, bp
+
+bool	pm_linenotempty()
+pointer	imgl2$t(), xt_fpval$t()
+
+begin
+	# If there are no bad pixels just get the image line and return.
+	if (fp == NULL)
+	    return (imgl2$t (im, line))
+
+	# Initialize
+	pm = FP_PM(fp)
+	nc = IM_LEN(im,1)
+	nl = IM_LEN(im,2)
+	ncols = FP_NCOLS(fp)
+	call amovkl (long(1), v, IM_MAXDIM)
+	v[2] = line
+
+	# If there might be column interpolation initialize value arrays.
+	if (ncols > 0 && FP_PV1(fp) == NULL) {
+	    FP_PIXTYPE(fp) = TY_PIXEL
+	    call malloc (FP_PV1(fp), ncols, FP_PIXTYPE(fp))
+	    call malloc (FP_PV2(fp), ncols, FP_PIXTYPE(fp))
+	    indef = INDEF
+	    call amovk$t (indef, Mem$t[FP_V1(fp,1)], ncols)
+	    call amovk$t (indef, Mem$t[FP_V2(fp,1)], ncols)
+	}
+
+	# If there are no bad pixels in the line and the line contains
+	# no column interpolation endpoints return the data directly.
+	# Otherwise get the line and fill in any endpoints that may
+	# be used later.
+
+	if (!pm_linenotempty (pm, v)) {
+	    if (line < FP_LMIN(fp) || line > FP_LMAX(fp))
+		return (imgl2$t (im, line))
+	    else
+		return (xt_fpval$t (fp, im, line))
+	}
+
+	# Get the pixel mask.
+	call malloc (bp, nc, TY_SHORT)
+	call pmglps (pm, v, Mems[bp], 0, nc, PIX_SRC)
+	bp = bp - 1
+
+	# Check if any column interpolation endpoints are needed and
+	# set them.  Set any other endpoints on the same lines at
+	# the same time.
+
+	if (line >= FP_LMIN(fp) && line < FP_LMAX(fp)) {
+	    j = 1
+	    do i = col1, col2 {
+		if (Mems[bp+i] == FP_CVAL(fp)) {
+		    for (; j<=ncols && FP_COL(fp,j)!=i; j=j+1)
+			;
+		    for (; j<=ncols && FP_COL(fp,j)==i; j=j+1) {
+			if (line>FP_L1(fp,j) && line<FP_L2(fp,j)) {
+			    if (IS_INDEF(Mem$t[FP_V1(fp,j)]))
+				data = xt_fpval$t (fp, im, FP_L1(fp,j))
+			    if (IS_INDEF(Mem$t[FP_V2(fp,j)]))
+				data = xt_fpval$t (fp, im, FP_L2(fp,j))
+			}
+		    }
+		}
+	    }
+	}
+
+	# Fix pixels by column or line interpolation.
+	if (FP_DATA(fp) == NULL) {
+	    FP_PIXTYPE(fp) = TY_PIXEL
+	    call malloc (FP_DATA(fp), nc, FP_PIXTYPE(fp))
+	}
+	data = FP_DATA(fp)
+	call amov$t (Mem$t[xt_fpval$t(fp,im,line)], Mem$t[data], nc)
+	j = 1
+	for (c1=col1; c1<=col2 && Mems[bp+c1]==0; c1=c1+1)
+	    ;
+	while (c1 <= col2) {
+	    c1 = c1 - 1
+	    for (c2=c1+2; c2<=col2 && Mems[bp+c2]!=0; c2=c2+1)
+		;
+	    a = INDEF
+	    do i = c1+1, c2-1 {
+		if (Mems[bp+i] == FP_LVAL(fp)) {
+		    if (IS_INDEF(a)) {
+			if (c1 < col1 && c2 > col2) {
+			    c1 = c2 + 1
+			    next
+			}
+			if (c1 >= col1)
+			    a = Mem$t[data+c1-1]
+			else
+			    a = Mem$t[data+c2-1]
+			if (c2 <= col2)
+			    b = (Mem$t[data+c2-1] - a) / (c2 - c1)
+			else
+			    b = 0.
+		    }
+		    val = a + b * (i - c1)
+		    if (fd != NULL) {
+			call fprintf (fd, "%4d %4d %8g %8g")
+			    call pargi (i)
+			    call pargi (line)
+			    call parg$t (Mem$t[data+i-1])
+			    $if (datatype == silr)
+			    call pargr (val)
+			    $else
+			    call parg$t (val)
+			    $endif
+			if (c1 >= col1) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (c1)
+			    call pargi (line)
+			}
+			if (c2 <= col2) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (c2)
+			    call pargi (line)
+			}
+			call fprintf (fd, "\n")
+		    }
+		} else {
+		    for (; j<=ncols && FP_COL(fp,j)!=i; j=j+1)
+			;
+		    for (; j<=ncols && FP_COL(fp,j)==i; j=j+1) {
+			l1 = FP_L1(fp,j)
+			l2 = FP_L2(fp,j)
+			if (l1 < line1 && l2 > line2)
+			    next
+			if (line > l1 && line < l2) {
+			    if (l1 >= line1)
+				c = Mem$t[FP_V1(fp,j)]
+			    else
+				c = Mem$t[FP_V2(fp,j)]
+			    if (l2 <= line2) {
+				d = (Mem$t[FP_V2(fp,j)] - c) / (l2 - l1)
+				val = c + d * (line - l1)
+			    } else
+				val = c
+			    l3 = l1
+			    l4 = l2
+			}
+		    }
+		    if (fd != NULL) {
+			call fprintf (fd, "%4d %4d %8g %8g")
+			    call pargi (i)
+			    call pargi (line)
+			    call parg$t (Mem$t[data+i-1])
+			    $if (datatype == silr)
+			    call pargr (val)
+			    $else
+			    call parg$t (val)
+			    $endif
+			if (l1 >= line1) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (i)
+			    call pargi (l3)
+			}
+			if (l2 <= line2) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (i)
+			    call pargi (l4)
+			}
+			call fprintf (fd, "\n")
+		    }
+		}
+		$if (datatype == sil)
+		Mem$t[data+i-1] = nint (val)
+		$else
+		Mem$t[data+i-1] = val
+		$endif
+	    }
+	    for (c1=c2+1; c1<=col2 && Mems[bp+c1]==0; c1=c1+1)
+		;
+	}
+
+	call mfree (bp, TY_SHORT)
+	return (data)
+end
+
+
+# XT_FPVAL -- Get data for the specified line and set the values for
+# all column interpolation endpoints which occur at that line.
+
+pointer procedure xt_fpval$t (fp, im, line)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+
+int	i
+pointer	data, imgl2$t()
+
+begin
+	# Set out of bounds values to 0.  These are not used but we need
+	# to cancel the INDEF values.
+	if (line < 1 || line > IM_LEN(im,2)) {
+	    do i = 1, FP_NCOLS(fp) {
+		if (line == FP_L1(fp,i))
+		    Mem$t[FP_V1(fp,i)] = 0.
+		else if (line == FP_L2(fp,i))
+		    Mem$t[FP_V2(fp,i)] = 0.
+	    }
+	    return (NULL)
+	}
+
+	data = imgl2$t (im, line)
+	do i = 1, FP_NCOLS(fp) {
+	    if (line == FP_L1(fp,i))
+		Mem$t[FP_V1(fp,i)] = Mem$t[data+FP_COL(fp,i)-1]
+	    else if (line == FP_L2(fp,i))
+		Mem$t[FP_V2(fp,i)] = Mem$t[data+FP_COL(fp,i)-1]
+	}
+
+	return (data)
+end
+
+$endfor
diff --git a/pkg/xtools/fixpix/xtfp.x b/pkg/xtools/fixpix/xtfp.x
new file mode 100644
index 00000000..774ffa12
--- /dev/null
+++ b/pkg/xtools/fixpix/xtfp.x
@@ -0,0 +1,1271 @@
+include	<imhdr.h>
+include	<pmset.h>
+include	"xtfixpix.h"
+
+
+
+
+# XT_FP -- Get the specified line of image data and replace bad pixels by
+# interpolation.
+
+pointer procedure xt_fps (fp, im, line, fd)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+int	fd			#I File descriptor for pixel list
+
+int	col1, col2		#I Section of interest
+int	line1, line2		#I Section of interest
+
+pointer	imgl2s(), xt_fpss()
+
+begin
+	# If there are no bad pixels just get the image line and return.
+	if (fp == NULL)
+	    return (imgl2s (im, line))
+
+	col1 = 1
+	col2 = IM_LEN(im,1)
+	line1 = 1
+	line2 = IM_LEN(im,2)
+
+	return (xt_fpss (fp, im, line, col1, col2, line1, line2, fd))
+end
+
+
+# XT_FXS -- Get the specified line of image data and replace bad pixels by
+# interpolation within a specified section.
+
+pointer procedure xt_fpss (fp, im, line, col1, col2, line1, line2, fd)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+int	fd			#I File descriptor for pixel list
+
+int	col1, col2		#I Section of interest
+int	line1, line2		#I Section of interest
+
+int	i, j, nc, nl, ncols, c1, c2, l1, l2, l3, l4
+long	v[IM_MAXDIM]
+real	a, b, c, d, val
+short	indef
+pointer	pm, data, bp
+
+bool	pm_linenotempty()
+pointer	imgl2s(), xt_fpvals()
+
+begin
+	# If there are no bad pixels just get the image line and return.
+	if (fp == NULL)
+	    return (imgl2s (im, line))
+
+	# Initialize
+	pm = FP_PM(fp)
+	nc = IM_LEN(im,1)
+	nl = IM_LEN(im,2)
+	ncols = FP_NCOLS(fp)
+	call amovkl (long(1), v, IM_MAXDIM)
+	v[2] = line
+
+	# If there might be column interpolation initialize value arrays.
+	if (ncols > 0 && FP_PV1(fp) == NULL) {
+	    FP_PIXTYPE(fp) = TY_SHORT
+	    call malloc (FP_PV1(fp), ncols, FP_PIXTYPE(fp))
+	    call malloc (FP_PV2(fp), ncols, FP_PIXTYPE(fp))
+	    indef = INDEFS
+	    call amovks (indef, Mems[FP_V1(fp,1)], ncols)
+	    call amovks (indef, Mems[FP_V2(fp,1)], ncols)
+	}
+
+	# If there are no bad pixels in the line and the line contains
+	# no column interpolation endpoints return the data directly.
+	# Otherwise get the line and fill in any endpoints that may
+	# be used later.
+
+	if (!pm_linenotempty (pm, v)) {
+	    if (line < FP_LMIN(fp) || line > FP_LMAX(fp))
+		return (imgl2s (im, line))
+	    else
+		return (xt_fpvals (fp, im, line))
+	}
+
+	# Get the pixel mask.
+	call malloc (bp, nc, TY_SHORT)
+	call pmglps (pm, v, Mems[bp], 0, nc, PIX_SRC)
+	bp = bp - 1
+
+	# Check if any column interpolation endpoints are needed and
+	# set them.  Set any other endpoints on the same lines at
+	# the same time.
+
+	if (line >= FP_LMIN(fp) && line < FP_LMAX(fp)) {
+	    j = 1
+	    do i = col1, col2 {
+		if (Mems[bp+i] == FP_CVAL(fp)) {
+		    for (; j<=ncols && FP_COL(fp,j)!=i; j=j+1)
+			;
+		    for (; j<=ncols && FP_COL(fp,j)==i; j=j+1) {
+			if (line>FP_L1(fp,j) && line<FP_L2(fp,j)) {
+			    if (IS_INDEFS(Mems[FP_V1(fp,j)]))
+				data = xt_fpvals (fp, im, FP_L1(fp,j))
+			    if (IS_INDEFS(Mems[FP_V2(fp,j)]))
+				data = xt_fpvals (fp, im, FP_L2(fp,j))
+			}
+		    }
+		}
+	    }
+	}
+
+	# Fix pixels by column or line interpolation.
+	if (FP_DATA(fp) == NULL) {
+	    FP_PIXTYPE(fp) = TY_SHORT
+	    call malloc (FP_DATA(fp), nc, FP_PIXTYPE(fp))
+	}
+	data = FP_DATA(fp)
+	call amovs (Mems[xt_fpvals(fp,im,line)], Mems[data], nc)
+	j = 1
+	for (c1=col1; c1<=col2 && Mems[bp+c1]==0; c1=c1+1)
+	    ;
+	while (c1 <= col2) {
+	    c1 = c1 - 1
+	    for (c2=c1+2; c2<=col2 && Mems[bp+c2]!=0; c2=c2+1)
+		;
+	    a = INDEFS
+	    do i = c1+1, c2-1 {
+		if (Mems[bp+i] == FP_LVAL(fp)) {
+		    if (IS_INDEFS(a)) {
+			if (c1 < col1 && c2 > col2) {
+			    c1 = c2 + 1
+			    next
+			}
+			if (c1 >= col1)
+			    a = Mems[data+c1-1]
+			else
+			    a = Mems[data+c2-1]
+			if (c2 <= col2)
+			    b = (Mems[data+c2-1] - a) / (c2 - c1)
+			else
+			    b = 0.
+		    }
+		    val = a + b * (i - c1)
+		    if (fd != NULL) {
+			call fprintf (fd, "%4d %4d %8g %8g")
+			    call pargi (i)
+			    call pargi (line)
+			    call pargs (Mems[data+i-1])
+			    call pargr (val)
+			if (c1 >= col1) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (c1)
+			    call pargi (line)
+			}
+			if (c2 <= col2) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (c2)
+			    call pargi (line)
+			}
+			call fprintf (fd, "\n")
+		    }
+		} else {
+		    for (; j<=ncols && FP_COL(fp,j)!=i; j=j+1)
+			;
+		    for (; j<=ncols && FP_COL(fp,j)==i; j=j+1) {
+			l1 = FP_L1(fp,j)
+			l2 = FP_L2(fp,j)
+			if (l1 < line1 && l2 > line2)
+			    next
+			if (line > l1 && line < l2) {
+			    if (l1 >= line1)
+				c = Mems[FP_V1(fp,j)]
+			    else
+				c = Mems[FP_V2(fp,j)]
+			    if (l2 <= line2) {
+				d = (Mems[FP_V2(fp,j)] - c) / (l2 - l1)
+				val = c + d * (line - l1)
+			    } else
+				val = c
+			    l3 = l1
+			    l4 = l2
+			}
+		    }
+		    if (fd != NULL) {
+			call fprintf (fd, "%4d %4d %8g %8g")
+			    call pargi (i)
+			    call pargi (line)
+			    call pargs (Mems[data+i-1])
+			    call pargr (val)
+			if (l1 >= line1) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (i)
+			    call pargi (l3)
+			}
+			if (l2 <= line2) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (i)
+			    call pargi (l4)
+			}
+			call fprintf (fd, "\n")
+		    }
+		}
+		Mems[data+i-1] = nint (val)
+	    }
+	    for (c1=c2+1; c1<=col2 && Mems[bp+c1]==0; c1=c1+1)
+		;
+	}
+
+	call mfree (bp, TY_SHORT)
+	return (data)
+end
+
+
+# XT_FPVAL -- Get data for the specified line and set the values for
+# all column interpolation endpoints which occur at that line.
+
+pointer procedure xt_fpvals (fp, im, line)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+
+int	i
+pointer	data, imgl2s()
+
+begin
+	# Set out of bounds values to 0.  These are not used but we need
+	# to cancel the INDEF values.
+	if (line < 1 || line > IM_LEN(im,2)) {
+	    do i = 1, FP_NCOLS(fp) {
+		if (line == FP_L1(fp,i))
+		    Mems[FP_V1(fp,i)] = 0.
+		else if (line == FP_L2(fp,i))
+		    Mems[FP_V2(fp,i)] = 0.
+	    }
+	    return (NULL)
+	}
+
+	data = imgl2s (im, line)
+	do i = 1, FP_NCOLS(fp) {
+	    if (line == FP_L1(fp,i))
+		Mems[FP_V1(fp,i)] = Mems[data+FP_COL(fp,i)-1]
+	    else if (line == FP_L2(fp,i))
+		Mems[FP_V2(fp,i)] = Mems[data+FP_COL(fp,i)-1]
+	}
+
+	return (data)
+end
+
+
+
+# XT_FP -- Get the specified line of image data and replace bad pixels by
+# interpolation.
+
+pointer procedure xt_fpi (fp, im, line, fd)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+int	fd			#I File descriptor for pixel list
+
+int	col1, col2		#I Section of interest
+int	line1, line2		#I Section of interest
+
+pointer	imgl2i(), xt_fpsi()
+
+begin
+	# If there are no bad pixels just get the image line and return.
+	if (fp == NULL)
+	    return (imgl2i (im, line))
+
+	col1 = 1
+	col2 = IM_LEN(im,1)
+	line1 = 1
+	line2 = IM_LEN(im,2)
+
+	return (xt_fpsi (fp, im, line, col1, col2, line1, line2, fd))
+end
+
+
+# XT_FXS -- Get the specified line of image data and replace bad pixels by
+# interpolation within a specified section.
+
+pointer procedure xt_fpsi (fp, im, line, col1, col2, line1, line2, fd)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+int	fd			#I File descriptor for pixel list
+
+int	col1, col2		#I Section of interest
+int	line1, line2		#I Section of interest
+
+int	i, j, nc, nl, ncols, c1, c2, l1, l2, l3, l4
+long	v[IM_MAXDIM]
+real	a, b, c, d, val
+int	indef
+pointer	pm, data, bp
+
+bool	pm_linenotempty()
+pointer	imgl2i(), xt_fpvali()
+
+begin
+	# If there are no bad pixels just get the image line and return.
+	if (fp == NULL)
+	    return (imgl2i (im, line))
+
+	# Initialize
+	pm = FP_PM(fp)
+	nc = IM_LEN(im,1)
+	nl = IM_LEN(im,2)
+	ncols = FP_NCOLS(fp)
+	call amovkl (long(1), v, IM_MAXDIM)
+	v[2] = line
+
+	# If there might be column interpolation initialize value arrays.
+	if (ncols > 0 && FP_PV1(fp) == NULL) {
+	    FP_PIXTYPE(fp) = TY_INT
+	    call malloc (FP_PV1(fp), ncols, FP_PIXTYPE(fp))
+	    call malloc (FP_PV2(fp), ncols, FP_PIXTYPE(fp))
+	    indef = INDEFI
+	    call amovki (indef, Memi[FP_V1(fp,1)], ncols)
+	    call amovki (indef, Memi[FP_V2(fp,1)], ncols)
+	}
+
+	# If there are no bad pixels in the line and the line contains
+	# no column interpolation endpoints return the data directly.
+	# Otherwise get the line and fill in any endpoints that may
+	# be used later.
+
+	if (!pm_linenotempty (pm, v)) {
+	    if (line < FP_LMIN(fp) || line > FP_LMAX(fp))
+		return (imgl2i (im, line))
+	    else
+		return (xt_fpvali (fp, im, line))
+	}
+
+	# Get the pixel mask.
+	call malloc (bp, nc, TY_SHORT)
+	call pmglps (pm, v, Mems[bp], 0, nc, PIX_SRC)
+	bp = bp - 1
+
+	# Check if any column interpolation endpoints are needed and
+	# set them.  Set any other endpoints on the same lines at
+	# the same time.
+
+	if (line >= FP_LMIN(fp) && line < FP_LMAX(fp)) {
+	    j = 1
+	    do i = col1, col2 {
+		if (Mems[bp+i] == FP_CVAL(fp)) {
+		    for (; j<=ncols && FP_COL(fp,j)!=i; j=j+1)
+			;
+		    for (; j<=ncols && FP_COL(fp,j)==i; j=j+1) {
+			if (line>FP_L1(fp,j) && line<FP_L2(fp,j)) {
+			    if (IS_INDEFI(Memi[FP_V1(fp,j)]))
+				data = xt_fpvali (fp, im, FP_L1(fp,j))
+			    if (IS_INDEFI(Memi[FP_V2(fp,j)]))
+				data = xt_fpvali (fp, im, FP_L2(fp,j))
+			}
+		    }
+		}
+	    }
+	}
+
+	# Fix pixels by column or line interpolation.
+	if (FP_DATA(fp) == NULL) {
+	    FP_PIXTYPE(fp) = TY_INT
+	    call malloc (FP_DATA(fp), nc, FP_PIXTYPE(fp))
+	}
+	data = FP_DATA(fp)
+	call amovi (Memi[xt_fpvali(fp,im,line)], Memi[data], nc)
+	j = 1
+	for (c1=col1; c1<=col2 && Mems[bp+c1]==0; c1=c1+1)
+	    ;
+	while (c1 <= col2) {
+	    c1 = c1 - 1
+	    for (c2=c1+2; c2<=col2 && Mems[bp+c2]!=0; c2=c2+1)
+		;
+	    a = INDEFI
+	    do i = c1+1, c2-1 {
+		if (Mems[bp+i] == FP_LVAL(fp)) {
+		    if (IS_INDEFI(a)) {
+			if (c1 < col1 && c2 > col2) {
+			    c1 = c2 + 1
+			    next
+			}
+			if (c1 >= col1)
+			    a = Memi[data+c1-1]
+			else
+			    a = Memi[data+c2-1]
+			if (c2 <= col2)
+			    b = (Memi[data+c2-1] - a) / (c2 - c1)
+			else
+			    b = 0.
+		    }
+		    val = a + b * (i - c1)
+		    if (fd != NULL) {
+			call fprintf (fd, "%4d %4d %8g %8g")
+			    call pargi (i)
+			    call pargi (line)
+			    call pargi (Memi[data+i-1])
+			    call pargr (val)
+			if (c1 >= col1) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (c1)
+			    call pargi (line)
+			}
+			if (c2 <= col2) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (c2)
+			    call pargi (line)
+			}
+			call fprintf (fd, "\n")
+		    }
+		} else {
+		    for (; j<=ncols && FP_COL(fp,j)!=i; j=j+1)
+			;
+		    for (; j<=ncols && FP_COL(fp,j)==i; j=j+1) {
+			l1 = FP_L1(fp,j)
+			l2 = FP_L2(fp,j)
+			if (l1 < line1 && l2 > line2)
+			    next
+			if (line > l1 && line < l2) {
+			    if (l1 >= line1)
+				c = Memi[FP_V1(fp,j)]
+			    else
+				c = Memi[FP_V2(fp,j)]
+			    if (l2 <= line2) {
+				d = (Memi[FP_V2(fp,j)] - c) / (l2 - l1)
+				val = c + d * (line - l1)
+			    } else
+				val = c
+			    l3 = l1
+			    l4 = l2
+			}
+		    }
+		    if (fd != NULL) {
+			call fprintf (fd, "%4d %4d %8g %8g")
+			    call pargi (i)
+			    call pargi (line)
+			    call pargi (Memi[data+i-1])
+			    call pargr (val)
+			if (l1 >= line1) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (i)
+			    call pargi (l3)
+			}
+			if (l2 <= line2) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (i)
+			    call pargi (l4)
+			}
+			call fprintf (fd, "\n")
+		    }
+		}
+		Memi[data+i-1] = nint (val)
+	    }
+	    for (c1=c2+1; c1<=col2 && Mems[bp+c1]==0; c1=c1+1)
+		;
+	}
+
+	call mfree (bp, TY_SHORT)
+	return (data)
+end
+
+
+# XT_FPVAL -- Get data for the specified line and set the values for
+# all column interpolation endpoints which occur at that line.
+
+pointer procedure xt_fpvali (fp, im, line)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+
+int	i
+pointer	data, imgl2i()
+
+begin
+	# Set out of bounds values to 0.  These are not used but we need
+	# to cancel the INDEF values.
+	if (line < 1 || line > IM_LEN(im,2)) {
+	    do i = 1, FP_NCOLS(fp) {
+		if (line == FP_L1(fp,i))
+		    Memi[FP_V1(fp,i)] = 0.
+		else if (line == FP_L2(fp,i))
+		    Memi[FP_V2(fp,i)] = 0.
+	    }
+	    return (NULL)
+	}
+
+	data = imgl2i (im, line)
+	do i = 1, FP_NCOLS(fp) {
+	    if (line == FP_L1(fp,i))
+		Memi[FP_V1(fp,i)] = Memi[data+FP_COL(fp,i)-1]
+	    else if (line == FP_L2(fp,i))
+		Memi[FP_V2(fp,i)] = Memi[data+FP_COL(fp,i)-1]
+	}
+
+	return (data)
+end
+
+
+
+# XT_FP -- Get the specified line of image data and replace bad pixels by
+# interpolation.
+
+pointer procedure xt_fpl (fp, im, line, fd)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+int	fd			#I File descriptor for pixel list
+
+int	col1, col2		#I Section of interest
+int	line1, line2		#I Section of interest
+
+pointer	imgl2l(), xt_fpsl()
+
+begin
+	# If there are no bad pixels just get the image line and return.
+	if (fp == NULL)
+	    return (imgl2l (im, line))
+
+	col1 = 1
+	col2 = IM_LEN(im,1)
+	line1 = 1
+	line2 = IM_LEN(im,2)
+
+	return (xt_fpsl (fp, im, line, col1, col2, line1, line2, fd))
+end
+
+
+# XT_FXS -- Get the specified line of image data and replace bad pixels by
+# interpolation within a specified section.
+
+pointer procedure xt_fpsl (fp, im, line, col1, col2, line1, line2, fd)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+int	fd			#I File descriptor for pixel list
+
+int	col1, col2		#I Section of interest
+int	line1, line2		#I Section of interest
+
+int	i, j, nc, nl, ncols, c1, c2, l1, l2, l3, l4
+long	v[IM_MAXDIM]
+real	a, b, c, d, val
+long	indef
+pointer	pm, data, bp
+
+bool	pm_linenotempty()
+pointer	imgl2l(), xt_fpvall()
+
+begin
+	# If there are no bad pixels just get the image line and return.
+	if (fp == NULL)
+	    return (imgl2l (im, line))
+
+	# Initialize
+	pm = FP_PM(fp)
+	nc = IM_LEN(im,1)
+	nl = IM_LEN(im,2)
+	ncols = FP_NCOLS(fp)
+	call amovkl (long(1), v, IM_MAXDIM)
+	v[2] = line
+
+	# If there might be column interpolation initialize value arrays.
+	if (ncols > 0 && FP_PV1(fp) == NULL) {
+	    FP_PIXTYPE(fp) = TY_LONG
+	    call malloc (FP_PV1(fp), ncols, FP_PIXTYPE(fp))
+	    call malloc (FP_PV2(fp), ncols, FP_PIXTYPE(fp))
+	    indef = INDEFL
+	    call amovkl (indef, Meml[FP_V1(fp,1)], ncols)
+	    call amovkl (indef, Meml[FP_V2(fp,1)], ncols)
+	}
+
+	# If there are no bad pixels in the line and the line contains
+	# no column interpolation endpoints return the data directly.
+	# Otherwise get the line and fill in any endpoints that may
+	# be used later.
+
+	if (!pm_linenotempty (pm, v)) {
+	    if (line < FP_LMIN(fp) || line > FP_LMAX(fp))
+		return (imgl2l (im, line))
+	    else
+		return (xt_fpvall (fp, im, line))
+	}
+
+	# Get the pixel mask.
+	call malloc (bp, nc, TY_SHORT)
+	call pmglps (pm, v, Mems[bp], 0, nc, PIX_SRC)
+	bp = bp - 1
+
+	# Check if any column interpolation endpoints are needed and
+	# set them.  Set any other endpoints on the same lines at
+	# the same time.
+
+	if (line >= FP_LMIN(fp) && line < FP_LMAX(fp)) {
+	    j = 1
+	    do i = col1, col2 {
+		if (Mems[bp+i] == FP_CVAL(fp)) {
+		    for (; j<=ncols && FP_COL(fp,j)!=i; j=j+1)
+			;
+		    for (; j<=ncols && FP_COL(fp,j)==i; j=j+1) {
+			if (line>FP_L1(fp,j) && line<FP_L2(fp,j)) {
+			    if (IS_INDEFL(Meml[FP_V1(fp,j)]))
+				data = xt_fpvall (fp, im, FP_L1(fp,j))
+			    if (IS_INDEFL(Meml[FP_V2(fp,j)]))
+				data = xt_fpvall (fp, im, FP_L2(fp,j))
+			}
+		    }
+		}
+	    }
+	}
+
+	# Fix pixels by column or line interpolation.
+	if (FP_DATA(fp) == NULL) {
+	    FP_PIXTYPE(fp) = TY_LONG
+	    call malloc (FP_DATA(fp), nc, FP_PIXTYPE(fp))
+	}
+	data = FP_DATA(fp)
+	call amovl (Meml[xt_fpvall(fp,im,line)], Meml[data], nc)
+	j = 1
+	for (c1=col1; c1<=col2 && Mems[bp+c1]==0; c1=c1+1)
+	    ;
+	while (c1 <= col2) {
+	    c1 = c1 - 1
+	    for (c2=c1+2; c2<=col2 && Mems[bp+c2]!=0; c2=c2+1)
+		;
+	    a = INDEFL
+	    do i = c1+1, c2-1 {
+		if (Mems[bp+i] == FP_LVAL(fp)) {
+		    if (IS_INDEFL(a)) {
+			if (c1 < col1 && c2 > col2) {
+			    c1 = c2 + 1
+			    next
+			}
+			if (c1 >= col1)
+			    a = Meml[data+c1-1]
+			else
+			    a = Meml[data+c2-1]
+			if (c2 <= col2)
+			    b = (Meml[data+c2-1] - a) / (c2 - c1)
+			else
+			    b = 0.
+		    }
+		    val = a + b * (i - c1)
+		    if (fd != NULL) {
+			call fprintf (fd, "%4d %4d %8g %8g")
+			    call pargi (i)
+			    call pargi (line)
+			    call pargl (Meml[data+i-1])
+			    call pargr (val)
+			if (c1 >= col1) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (c1)
+			    call pargi (line)
+			}
+			if (c2 <= col2) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (c2)
+			    call pargi (line)
+			}
+			call fprintf (fd, "\n")
+		    }
+		} else {
+		    for (; j<=ncols && FP_COL(fp,j)!=i; j=j+1)
+			;
+		    for (; j<=ncols && FP_COL(fp,j)==i; j=j+1) {
+			l1 = FP_L1(fp,j)
+			l2 = FP_L2(fp,j)
+			if (l1 < line1 && l2 > line2)
+			    next
+			if (line > l1 && line < l2) {
+			    if (l1 >= line1)
+				c = Meml[FP_V1(fp,j)]
+			    else
+				c = Meml[FP_V2(fp,j)]
+			    if (l2 <= line2) {
+				d = (Meml[FP_V2(fp,j)] - c) / (l2 - l1)
+				val = c + d * (line - l1)
+			    } else
+				val = c
+			    l3 = l1
+			    l4 = l2
+			}
+		    }
+		    if (fd != NULL) {
+			call fprintf (fd, "%4d %4d %8g %8g")
+			    call pargi (i)
+			    call pargi (line)
+			    call pargl (Meml[data+i-1])
+			    call pargr (val)
+			if (l1 >= line1) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (i)
+			    call pargi (l3)
+			}
+			if (l2 <= line2) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (i)
+			    call pargi (l4)
+			}
+			call fprintf (fd, "\n")
+		    }
+		}
+		Meml[data+i-1] = nint (val)
+	    }
+	    for (c1=c2+1; c1<=col2 && Mems[bp+c1]==0; c1=c1+1)
+		;
+	}
+
+	call mfree (bp, TY_SHORT)
+	return (data)
+end
+
+
+# XT_FPVAL -- Get data for the specified line and set the values for
+# all column interpolation endpoints which occur at that line.
+
+pointer procedure xt_fpvall (fp, im, line)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+
+int	i
+pointer	data, imgl2l()
+
+begin
+	# Set out of bounds values to 0.  These are not used but we need
+	# to cancel the INDEF values.
+	if (line < 1 || line > IM_LEN(im,2)) {
+	    do i = 1, FP_NCOLS(fp) {
+		if (line == FP_L1(fp,i))
+		    Meml[FP_V1(fp,i)] = 0.
+		else if (line == FP_L2(fp,i))
+		    Meml[FP_V2(fp,i)] = 0.
+	    }
+	    return (NULL)
+	}
+
+	data = imgl2l (im, line)
+	do i = 1, FP_NCOLS(fp) {
+	    if (line == FP_L1(fp,i))
+		Meml[FP_V1(fp,i)] = Meml[data+FP_COL(fp,i)-1]
+	    else if (line == FP_L2(fp,i))
+		Meml[FP_V2(fp,i)] = Meml[data+FP_COL(fp,i)-1]
+	}
+
+	return (data)
+end
+
+
+
+# XT_FP -- Get the specified line of image data and replace bad pixels by
+# interpolation.
+
+pointer procedure xt_fpr (fp, im, line, fd)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+int	fd			#I File descriptor for pixel list
+
+int	col1, col2		#I Section of interest
+int	line1, line2		#I Section of interest
+
+pointer	imgl2r(), xt_fpsr()
+
+begin
+	# If there are no bad pixels just get the image line and return.
+	if (fp == NULL)
+	    return (imgl2r (im, line))
+
+	col1 = 1
+	col2 = IM_LEN(im,1)
+	line1 = 1
+	line2 = IM_LEN(im,2)
+
+	return (xt_fpsr (fp, im, line, col1, col2, line1, line2, fd))
+end
+
+
+# XT_FXS -- Get the specified line of image data and replace bad pixels by
+# interpolation within a specified section.
+
+pointer procedure xt_fpsr (fp, im, line, col1, col2, line1, line2, fd)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+int	fd			#I File descriptor for pixel list
+
+int	col1, col2		#I Section of interest
+int	line1, line2		#I Section of interest
+
+int	i, j, nc, nl, ncols, c1, c2, l1, l2, l3, l4
+long	v[IM_MAXDIM]
+real	a, b, c, d, val
+real	indef
+pointer	pm, data, bp
+
+bool	pm_linenotempty()
+pointer	imgl2r(), xt_fpvalr()
+
+begin
+	# If there are no bad pixels just get the image line and return.
+	if (fp == NULL)
+	    return (imgl2r (im, line))
+
+	# Initialize
+	pm = FP_PM(fp)
+	nc = IM_LEN(im,1)
+	nl = IM_LEN(im,2)
+	ncols = FP_NCOLS(fp)
+	call amovkl (long(1), v, IM_MAXDIM)
+	v[2] = line
+
+	# If there might be column interpolation initialize value arrays.
+	if (ncols > 0 && FP_PV1(fp) == NULL) {
+	    FP_PIXTYPE(fp) = TY_REAL
+	    call malloc (FP_PV1(fp), ncols, FP_PIXTYPE(fp))
+	    call malloc (FP_PV2(fp), ncols, FP_PIXTYPE(fp))
+	    indef = INDEFR
+	    call amovkr (indef, Memr[FP_V1(fp,1)], ncols)
+	    call amovkr (indef, Memr[FP_V2(fp,1)], ncols)
+	}
+
+	# If there are no bad pixels in the line and the line contains
+	# no column interpolation endpoints return the data directly.
+	# Otherwise get the line and fill in any endpoints that may
+	# be used later.
+
+	if (!pm_linenotempty (pm, v)) {
+	    if (line < FP_LMIN(fp) || line > FP_LMAX(fp))
+		return (imgl2r (im, line))
+	    else
+		return (xt_fpvalr (fp, im, line))
+	}
+
+	# Get the pixel mask.
+	call malloc (bp, nc, TY_SHORT)
+	call pmglps (pm, v, Mems[bp], 0, nc, PIX_SRC)
+	bp = bp - 1
+
+	# Check if any column interpolation endpoints are needed and
+	# set them.  Set any other endpoints on the same lines at
+	# the same time.
+
+	if (line >= FP_LMIN(fp) && line < FP_LMAX(fp)) {
+	    j = 1
+	    do i = col1, col2 {
+		if (Mems[bp+i] == FP_CVAL(fp)) {
+		    for (; j<=ncols && FP_COL(fp,j)!=i; j=j+1)
+			;
+		    for (; j<=ncols && FP_COL(fp,j)==i; j=j+1) {
+			if (line>FP_L1(fp,j) && line<FP_L2(fp,j)) {
+			    if (IS_INDEFR(Memr[FP_V1(fp,j)]))
+				data = xt_fpvalr (fp, im, FP_L1(fp,j))
+			    if (IS_INDEFR(Memr[FP_V2(fp,j)]))
+				data = xt_fpvalr (fp, im, FP_L2(fp,j))
+			}
+		    }
+		}
+	    }
+	}
+
+	# Fix pixels by column or line interpolation.
+	if (FP_DATA(fp) == NULL) {
+	    FP_PIXTYPE(fp) = TY_REAL
+	    call malloc (FP_DATA(fp), nc, FP_PIXTYPE(fp))
+	}
+	data = FP_DATA(fp)
+	call amovr (Memr[xt_fpvalr(fp,im,line)], Memr[data], nc)
+	j = 1
+	for (c1=col1; c1<=col2 && Mems[bp+c1]==0; c1=c1+1)
+	    ;
+	while (c1 <= col2) {
+	    c1 = c1 - 1
+	    for (c2=c1+2; c2<=col2 && Mems[bp+c2]!=0; c2=c2+1)
+		;
+	    a = INDEFR
+	    do i = c1+1, c2-1 {
+		if (Mems[bp+i] == FP_LVAL(fp)) {
+		    if (IS_INDEFR(a)) {
+			if (c1 < col1 && c2 > col2) {
+			    c1 = c2 + 1
+			    next
+			}
+			if (c1 >= col1)
+			    a = Memr[data+c1-1]
+			else
+			    a = Memr[data+c2-1]
+			if (c2 <= col2)
+			    b = (Memr[data+c2-1] - a) / (c2 - c1)
+			else
+			    b = 0.
+		    }
+		    val = a + b * (i - c1)
+		    if (fd != NULL) {
+			call fprintf (fd, "%4d %4d %8g %8g")
+			    call pargi (i)
+			    call pargi (line)
+			    call pargr (Memr[data+i-1])
+			    call pargr (val)
+			if (c1 >= col1) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (c1)
+			    call pargi (line)
+			}
+			if (c2 <= col2) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (c2)
+			    call pargi (line)
+			}
+			call fprintf (fd, "\n")
+		    }
+		} else {
+		    for (; j<=ncols && FP_COL(fp,j)!=i; j=j+1)
+			;
+		    for (; j<=ncols && FP_COL(fp,j)==i; j=j+1) {
+			l1 = FP_L1(fp,j)
+			l2 = FP_L2(fp,j)
+			if (l1 < line1 && l2 > line2)
+			    next
+			if (line > l1 && line < l2) {
+			    if (l1 >= line1)
+				c = Memr[FP_V1(fp,j)]
+			    else
+				c = Memr[FP_V2(fp,j)]
+			    if (l2 <= line2) {
+				d = (Memr[FP_V2(fp,j)] - c) / (l2 - l1)
+				val = c + d * (line - l1)
+			    } else
+				val = c
+			    l3 = l1
+			    l4 = l2
+			}
+		    }
+		    if (fd != NULL) {
+			call fprintf (fd, "%4d %4d %8g %8g")
+			    call pargi (i)
+			    call pargi (line)
+			    call pargr (Memr[data+i-1])
+			    call pargr (val)
+			if (l1 >= line1) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (i)
+			    call pargi (l3)
+			}
+			if (l2 <= line2) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (i)
+			    call pargi (l4)
+			}
+			call fprintf (fd, "\n")
+		    }
+		}
+		Memr[data+i-1] = val
+	    }
+	    for (c1=c2+1; c1<=col2 && Mems[bp+c1]==0; c1=c1+1)
+		;
+	}
+
+	call mfree (bp, TY_SHORT)
+	return (data)
+end
+
+
+# XT_FPVAL -- Get data for the specified line and set the values for
+# all column interpolation endpoints which occur at that line.
+
+pointer procedure xt_fpvalr (fp, im, line)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+
+int	i
+pointer	data, imgl2r()
+
+begin
+	# Set out of bounds values to 0.  These are not used but we need
+	# to cancel the INDEF values.
+	if (line < 1 || line > IM_LEN(im,2)) {
+	    do i = 1, FP_NCOLS(fp) {
+		if (line == FP_L1(fp,i))
+		    Memr[FP_V1(fp,i)] = 0.
+		else if (line == FP_L2(fp,i))
+		    Memr[FP_V2(fp,i)] = 0.
+	    }
+	    return (NULL)
+	}
+
+	data = imgl2r (im, line)
+	do i = 1, FP_NCOLS(fp) {
+	    if (line == FP_L1(fp,i))
+		Memr[FP_V1(fp,i)] = Memr[data+FP_COL(fp,i)-1]
+	    else if (line == FP_L2(fp,i))
+		Memr[FP_V2(fp,i)] = Memr[data+FP_COL(fp,i)-1]
+	}
+
+	return (data)
+end
+
+
+
+# XT_FP -- Get the specified line of image data and replace bad pixels by
+# interpolation.
+
+pointer procedure xt_fpd (fp, im, line, fd)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+int	fd			#I File descriptor for pixel list
+
+int	col1, col2		#I Section of interest
+int	line1, line2		#I Section of interest
+
+pointer	imgl2d(), xt_fpsd()
+
+begin
+	# If there are no bad pixels just get the image line and return.
+	if (fp == NULL)
+	    return (imgl2d (im, line))
+
+	col1 = 1
+	col2 = IM_LEN(im,1)
+	line1 = 1
+	line2 = IM_LEN(im,2)
+
+	return (xt_fpsd (fp, im, line, col1, col2, line1, line2, fd))
+end
+
+
+# XT_FXS -- Get the specified line of image data and replace bad pixels by
+# interpolation within a specified section.
+
+pointer procedure xt_fpsd (fp, im, line, col1, col2, line1, line2, fd)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+int	fd			#I File descriptor for pixel list
+
+int	col1, col2		#I Section of interest
+int	line1, line2		#I Section of interest
+
+int	i, j, nc, nl, ncols, c1, c2, l1, l2, l3, l4
+long	v[IM_MAXDIM]
+double	a, b, c, d, val
+double	indef
+pointer	pm, data, bp
+
+bool	pm_linenotempty()
+pointer	imgl2d(), xt_fpvald()
+
+begin
+	# If there are no bad pixels just get the image line and return.
+	if (fp == NULL)
+	    return (imgl2d (im, line))
+
+	# Initialize
+	pm = FP_PM(fp)
+	nc = IM_LEN(im,1)
+	nl = IM_LEN(im,2)
+	ncols = FP_NCOLS(fp)
+	call amovkl (long(1), v, IM_MAXDIM)
+	v[2] = line
+
+	# If there might be column interpolation initialize value arrays.
+	if (ncols > 0 && FP_PV1(fp) == NULL) {
+	    FP_PIXTYPE(fp) = TY_DOUBLE
+	    call malloc (FP_PV1(fp), ncols, FP_PIXTYPE(fp))
+	    call malloc (FP_PV2(fp), ncols, FP_PIXTYPE(fp))
+	    indef = INDEFD
+	    call amovkd (indef, Memd[FP_V1(fp,1)], ncols)
+	    call amovkd (indef, Memd[FP_V2(fp,1)], ncols)
+	}
+
+	# If there are no bad pixels in the line and the line contains
+	# no column interpolation endpoints return the data directly.
+	# Otherwise get the line and fill in any endpoints that may
+	# be used later.
+
+	if (!pm_linenotempty (pm, v)) {
+	    if (line < FP_LMIN(fp) || line > FP_LMAX(fp))
+		return (imgl2d (im, line))
+	    else
+		return (xt_fpvald (fp, im, line))
+	}
+
+	# Get the pixel mask.
+	call malloc (bp, nc, TY_SHORT)
+	call pmglps (pm, v, Mems[bp], 0, nc, PIX_SRC)
+	bp = bp - 1
+
+	# Check if any column interpolation endpoints are needed and
+	# set them.  Set any other endpoints on the same lines at
+	# the same time.
+
+	if (line >= FP_LMIN(fp) && line < FP_LMAX(fp)) {
+	    j = 1
+	    do i = col1, col2 {
+		if (Mems[bp+i] == FP_CVAL(fp)) {
+		    for (; j<=ncols && FP_COL(fp,j)!=i; j=j+1)
+			;
+		    for (; j<=ncols && FP_COL(fp,j)==i; j=j+1) {
+			if (line>FP_L1(fp,j) && line<FP_L2(fp,j)) {
+			    if (IS_INDEFD(Memd[FP_V1(fp,j)]))
+				data = xt_fpvald (fp, im, FP_L1(fp,j))
+			    if (IS_INDEFD(Memd[FP_V2(fp,j)]))
+				data = xt_fpvald (fp, im, FP_L2(fp,j))
+			}
+		    }
+		}
+	    }
+	}
+
+	# Fix pixels by column or line interpolation.
+	if (FP_DATA(fp) == NULL) {
+	    FP_PIXTYPE(fp) = TY_DOUBLE
+	    call malloc (FP_DATA(fp), nc, FP_PIXTYPE(fp))
+	}
+	data = FP_DATA(fp)
+	call amovd (Memd[xt_fpvald(fp,im,line)], Memd[data], nc)
+	j = 1
+	for (c1=col1; c1<=col2 && Mems[bp+c1]==0; c1=c1+1)
+	    ;
+	while (c1 <= col2) {
+	    c1 = c1 - 1
+	    for (c2=c1+2; c2<=col2 && Mems[bp+c2]!=0; c2=c2+1)
+		;
+	    a = INDEFD
+	    do i = c1+1, c2-1 {
+		if (Mems[bp+i] == FP_LVAL(fp)) {
+		    if (IS_INDEFD(a)) {
+			if (c1 < col1 && c2 > col2) {
+			    c1 = c2 + 1
+			    next
+			}
+			if (c1 >= col1)
+			    a = Memd[data+c1-1]
+			else
+			    a = Memd[data+c2-1]
+			if (c2 <= col2)
+			    b = (Memd[data+c2-1] - a) / (c2 - c1)
+			else
+			    b = 0.
+		    }
+		    val = a + b * (i - c1)
+		    if (fd != NULL) {
+			call fprintf (fd, "%4d %4d %8g %8g")
+			    call pargi (i)
+			    call pargi (line)
+			    call pargd (Memd[data+i-1])
+			    call pargd (val)
+			if (c1 >= col1) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (c1)
+			    call pargi (line)
+			}
+			if (c2 <= col2) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (c2)
+			    call pargi (line)
+			}
+			call fprintf (fd, "\n")
+		    }
+		} else {
+		    for (; j<=ncols && FP_COL(fp,j)!=i; j=j+1)
+			;
+		    for (; j<=ncols && FP_COL(fp,j)==i; j=j+1) {
+			l1 = FP_L1(fp,j)
+			l2 = FP_L2(fp,j)
+			if (l1 < line1 && l2 > line2)
+			    next
+			if (line > l1 && line < l2) {
+			    if (l1 >= line1)
+				c = Memd[FP_V1(fp,j)]
+			    else
+				c = Memd[FP_V2(fp,j)]
+			    if (l2 <= line2) {
+				d = (Memd[FP_V2(fp,j)] - c) / (l2 - l1)
+				val = c + d * (line - l1)
+			    } else
+				val = c
+			    l3 = l1
+			    l4 = l2
+			}
+		    }
+		    if (fd != NULL) {
+			call fprintf (fd, "%4d %4d %8g %8g")
+			    call pargi (i)
+			    call pargi (line)
+			    call pargd (Memd[data+i-1])
+			    call pargd (val)
+			if (l1 >= line1) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (i)
+			    call pargi (l3)
+			}
+			if (l2 <= line2) {
+			    call fprintf (fd, " %4d %4d")
+			    call pargi (i)
+			    call pargi (l4)
+			}
+			call fprintf (fd, "\n")
+		    }
+		}
+		Memd[data+i-1] = val
+	    }
+	    for (c1=c2+1; c1<=col2 && Mems[bp+c1]==0; c1=c1+1)
+		;
+	}
+
+	call mfree (bp, TY_SHORT)
+	return (data)
+end
+
+
+# XT_FPVAL -- Get data for the specified line and set the values for
+# all column interpolation endpoints which occur at that line.
+
+pointer procedure xt_fpvald (fp, im, line)
+
+pointer	fp			#I FIXPIX pointer
+pointer	im			#I Image pointer
+int	line			#I Line
+
+int	i
+pointer	data, imgl2d()
+
+begin
+	# Set out of bounds values to 0.  These are not used but we need
+	# to cancel the INDEF values.
+	if (line < 1 || line > IM_LEN(im,2)) {
+	    do i = 1, FP_NCOLS(fp) {
+		if (line == FP_L1(fp,i))
+		    Memd[FP_V1(fp,i)] = 0.
+		else if (line == FP_L2(fp,i))
+		    Memd[FP_V2(fp,i)] = 0.
+	    }
+	    return (NULL)
+	}
+
+	data = imgl2d (im, line)
+	do i = 1, FP_NCOLS(fp) {
+	    if (line == FP_L1(fp,i))
+		Memd[FP_V1(fp,i)] = Memd[data+FP_COL(fp,i)-1]
+	    else if (line == FP_L2(fp,i))
+		Memd[FP_V2(fp,i)] = Memd[data+FP_COL(fp,i)-1]
+	}
+
+	return (data)
+end
+
+
diff --git a/pkg/xtools/fixpix/xtpmmap.x b/pkg/xtools/fixpix/xtpmmap.x
new file mode 100644
index 00000000..54bbf954
--- /dev/null
+++ b/pkg/xtools/fixpix/xtpmmap.x
@@ -0,0 +1,693 @@
+include	<mach.h>
+include	<ctype.h>
+include	<error.h>
+include	<imhdr.h>
+include	<imset.h>
+include	<pmset.h>
+include	<mwset.h>
+include	<syserr.h>
+
+
+# XT_PMMAP -- Open a pixel mask READ_ONLY.
+#
+# This routine maps multiple types of mask files and designations.
+# It matches the mask coordinates to the reference image based on the
+# physical coordinate system so the mask may be of a different size.
+# The mask name is returned so that the task has the name pointed to by "BPM".
+# A null filename is allowed and returns NULL.
+
+pointer procedure xt_pmmap (pmname, refim, mname, sz_mname)
+
+char	pmname[ARB]		#I Pixel mask name
+pointer	refim			#I Reference image pointer
+char	mname[ARB]		#O Expanded mask name
+int	sz_mname		#O Size of expanded mask name
+
+int	i, flag, nowhite()
+pointer	sp, fname, im, ref, xt_pmmap1()
+bool	streq()
+errchk	xt_pmmap1
+
+begin
+	call smark (sp)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+
+	im = NULL
+	i = nowhite (pmname, Memc[fname], SZ_FNAME)
+	if (Memc[fname] == '!') {
+	    iferr (call imgstr (refim, Memc[fname+1], Memc[fname], SZ_FNAME))
+		Memc[fname] = EOS
+	} else if (streq (Memc[fname], "BPM")) {
+	    iferr (call imgstr (refim, "BPM", Memc[fname], SZ_FNAME))
+		Memc[fname] = EOS
+	} else if (streq (Memc[fname], "^BPM")) {
+	    flag = INVERT_MASK
+	    iferr (call imgstr (refim, "BPM", Memc[fname+1], SZ_FNAME))
+		Memc[fname] = EOS
+	}
+
+	if (Memc[fname] == '^') {
+	    flag = INVERT_MASK
+	    call strcpy (Memc[fname+1], Memc[fname], SZ_FNAME)
+	} else
+	    flag = NO
+
+	if (streq (Memc[fname], "EMPTY"))
+	    ref = refim
+	else
+	    ref = NULL
+
+	if (Memc[fname] != EOS)
+	    im = xt_pmmap1 (Memc[fname], ref, refim, flag, YES)
+	call strcpy (Memc[fname], mname, sz_mname)
+
+	call sfree (sp)
+	return (im)
+end
+
+
+# XT_MAPPM -- Open a pixel mask READ_ONLY with/without matching.
+#
+# This routine maps multiple types of mask files and designations.
+# It may match the mask coordinates to the reference image based on the
+# physical coordinate system.  In either case the mask is matched to be
+# the same size.  The mask name is returned so that the task has the
+# name pointed to by "BPM".  A null filename is allowed and returns NULL.
+
+pointer procedure xt_mappm (pmname, refim, match, mname, sz_mname)
+
+char	pmname[ARB]		#I Pixel mask name
+pointer	refim			#I Reference image pointer
+int	match			#I Match by physical coordinates?
+char	mname[ARB]		#O Expanded mask name
+int	sz_mname		#O Size of expanded mask name
+
+int	i, flag, nowhite()
+pointer	sp, fname, im, ref, xt_pmmap1()
+bool	streq()
+errchk	xt_pmmap1
+
+begin
+	call smark (sp)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+
+	im = NULL
+	i = nowhite (pmname, Memc[fname], SZ_FNAME)
+	if (Memc[fname] == '!') {
+	    iferr (call imgstr (refim, Memc[fname+1], Memc[fname], SZ_FNAME))
+		Memc[fname] = EOS
+	} else if (streq (Memc[fname], "BPM")) {
+	    iferr (call imgstr (refim, "BPM", Memc[fname], SZ_FNAME))
+		Memc[fname] = EOS
+	} else if (streq (Memc[fname], "^BPM")) {
+	    flag = INVERT_MASK
+	    iferr (call imgstr (refim, "BPM", Memc[fname+1], SZ_FNAME))
+		Memc[fname] = EOS
+	}
+
+	if (Memc[fname] == '^') {
+	    flag = INVERT_MASK
+	    call strcpy (Memc[fname+1], Memc[fname], SZ_FNAME)
+	} else
+	    flag = NO
+
+	if (streq (Memc[fname], "EMPTY"))
+	    ref = refim
+	else
+	    ref = NULL
+
+	if (Memc[fname] != EOS)
+	    im = xt_pmmap1 (Memc[fname], ref, refim, flag, match)
+	call strcpy (Memc[fname], mname, sz_mname)
+
+	call sfree (sp)
+	return (im)
+end
+
+
+# XT_PMUNMAP -- Unmap a mask image.
+# Note that the imio pointer may be purely an internal pointer opened
+# with im_pmmapo so we need to free the pl pointer explicitly.
+
+procedure xt_pmunmap (im)
+
+pointer	im			#I IMIO pointer for mask
+
+pointer	pm
+int	imstati()
+
+begin
+	pm = imstati (im, IM_PMDES)
+	call pm_close (pm)
+	call imseti (im, IM_PMDES, NULL)
+	call imunmap (im)
+end
+
+
+# XT_PMMAP1 -- Open a pixel mask READ_ONLY.  The input mask may be
+# a pixel list image, a non-pixel list image, or a text file.
+# Return error if the pixel mask cannot be opened.  For pixel masks
+# or image masks match the WCS.
+
+pointer procedure xt_pmmap1 (pmname, ref, refim, flag, match)
+
+char	pmname[ARB]		#I Pixel mask name
+pointer	ref			#I Reference image for pixel mask
+pointer	refim			#I Reference image for image or text
+int	flag			#I Mask flag
+int	match			#I Match by physical coordinates?
+
+int	imstati(),  errcode()
+pointer	im, pm
+pointer	im_pmmap(), xt_pmimmap(), xt_pmtext(), xt_pmsection()
+bool	streq()
+errchk	xt_match
+
+begin
+	im = NULL
+
+	if (streq (pmname, "STDIN"))
+	    im = xt_pmtext (pmname, refim, flag)
+
+	else if (pmname[1] == '[')
+	    im = xt_pmsection (pmname, refim, flag)
+
+	else {
+	    ifnoerr (im = im_pmmap (pmname, READ_ONLY, ref)) {
+		call xt_match (im, refim, match)
+		if (flag == INVERT_MASK) {
+		    pm = imstati (im, IM_PMDES)
+		    call xt_pminvert (pm)
+		    call imseti (im, IM_PMDES, pm)
+		}
+	    } else {
+		switch (errcode()) {
+		case SYS_IKIOPEN, SYS_FOPNNEXFIL, SYS_PLBADSAVEF, SYS_FOPEN:
+		    ifnoerr (im = xt_pmimmap (pmname, refim, flag))
+			call xt_match (im, refim, match)
+		    else {
+			switch (errcode()) {
+			case SYS_IKIOPEN:
+			    im = xt_pmtext (pmname, refim, flag)
+			default:
+			    call erract (EA_ERROR)
+			}
+		    }
+		default:
+		    call erract (EA_ERROR)
+		}
+	    }
+	}
+
+	return (im)
+end
+
+
+# XT_PMIMMAP -- Open a pixel mask from a non-pixel list image.
+# Return error if the image cannot be opened.
+
+pointer procedure xt_pmimmap (pmname, refim, flag)
+
+char	pmname[ARB]		#I Image name
+pointer	refim			#I Reference image pointer
+int	flag			#I Mask flag
+
+int	i, ndim, npix, rop, val
+pointer	sp, v1, v2, im_in, im_out, pm, mw, data
+
+int	imstati(), imgnli()
+pointer immap(), pm_newmask(), im_pmmapo(), imgl1i(), mw_openim()
+errchk	immap, mw_openim, im_pmmapo
+
+begin
+	call smark (sp)
+	call salloc (v1, IM_MAXDIM, TY_LONG)
+	call salloc (v2, IM_MAXDIM, TY_LONG)
+
+	call amovkl (long(1), Meml[v1], IM_MAXDIM)
+	call amovkl (long(1), Meml[v2], IM_MAXDIM)
+
+	im_in = immap (pmname, READ_ONLY, 0)
+	pm = imstati (im_in, IM_PMDES)
+	if (pm != NULL)
+	    return (im_in)
+	pm = pm_newmask (im_in, 16)
+
+	ndim = IM_NDIM(im_in)
+	npix = IM_LEN(im_in,1)
+
+	if (flag == INVERT_MASK)
+	    rop = PIX_NOT(PIX_SRC)
+	else
+	    rop = PIX_SRC
+
+	while (imgnli (im_in, data, Meml[v1]) != EOF) {
+	    if (flag == INVERT_MASK) {
+		do i = 0, npix-1 {
+		    val = Memi[data+i]
+		    if (val <= 0)
+		       Memi[data+i] = 1
+		    else
+		       Memi[data+i] = 0
+		}
+	    } else {
+		do i = 0, npix-1 {
+		    val = Memi[data+i]
+		    if (val < 0)
+		       Memi[data+i] = 0
+		}
+	    }
+	    call pmplpi (pm, Meml[v2], Memi[data], 0, npix, rop)
+	    call amovl (Meml[v1], Meml[v2], ndim)
+	}
+
+	im_out = im_pmmapo (pm, im_in)
+	data = imgl1i (im_out)		# Force I/O to set header
+	mw = mw_openim (im_in)		# Set WCS
+	call mw_saveim (mw, im_out)
+	call mw_close (mw)
+
+	#call imunmap (im_in)
+	call xt_pmunmap (im_in)
+	call sfree (sp)
+	return (im_out)
+end
+
+
+# XT_PMTEXT -- Create a pixel mask from a text file of rectangles.
+# Return error if the file cannot be opened.
+# This routine only applies to the first 2D plane.
+
+pointer procedure xt_pmtext (pmname, refim, flag)
+
+char	pmname[ARB]		#I Image name
+pointer	refim			#I Reference image pointer
+int	flag			#I Mask flag
+
+int	fd, nc, nl, c1, c2, l1, l2, nc1, nl1, rop
+pointer	pm, im, mw, dummy
+
+int	open(), fscan(), nscan()
+pointer	pm_newmask(), im_pmmapo(), imgl1i(), mw_openim()
+errchk	open, im_pmmapo
+
+begin
+	fd = open (pmname, READ_ONLY, TEXT_FILE)
+	pm = pm_newmask (refim, 16)
+
+	nc = IM_LEN(refim,1)
+	nl = IM_LEN(refim,2)
+
+	if (flag == INVERT_MASK)
+	    call pl_box (pm, 1, 1, nc, nl, PIX_SET+PIX_VALUE(1))
+
+	while (fscan (fd) != EOF) {
+	    call gargi (c1)
+	    call gargi (c2)
+	    call gargi (l1)
+	    call gargi (l2)
+	    if (nscan() != 4) {
+		if (nscan() == 2) {
+		    l1 = c2
+		    c2 = c1
+		    l2 = l1
+		} else
+		    next
+	    }
+
+	    c1 = max (1, c1)
+	    c2 = min (nc, c2)
+	    l1 = max (1, l1)
+	    l2 = min (nl, l2)
+	    nc1 = c2 - c1 + 1
+	    nl1 = l2 - l1 + 1
+	    if (nc1 < 1 || nl1 < 1)
+		next
+
+	    # Select mask value based on shape of rectangle.
+	    if (flag == INVERT_MASK)
+		rop = PIX_CLR
+	    else if (nc1 <= nl1)
+		rop = PIX_SET+PIX_VALUE(2)
+	    else
+		rop = PIX_SET+PIX_VALUE(3)
+
+	    # Set mask rectangle.
+	    call pm_box (pm, c1, l1, c2, l2, rop)
+	}
+
+	call close (fd)
+	im = im_pmmapo (pm, refim)
+	dummy = imgl1i (im)		# Force I/O to set header
+	mw = mw_openim (refim)		# Set WCS
+	call mw_saveim (mw, im)
+	call mw_close (mw)
+
+	return (im)
+end
+
+
+# XT_PMSECTION -- Create a pixel mask from an image section.
+# This only applies the mask to the first plane of the image.
+
+pointer procedure xt_pmsection (section, refim, flag)
+
+char	section[ARB]		#I Image section
+pointer	refim			#I Reference image pointer
+int	flag			#I Mask flag
+
+int	i, j, ip, temp, a[2], b[2], c[2], rop, ctoi()
+pointer	pm, im, mw, dummy, pm_newmask(), im_pmmapo(), imgl1i(), mw_openim()
+errchk	im_pmmapo
+define  error_  99
+
+begin
+	# This is currently only for 1D and 2D images.
+	if (IM_NDIM(refim) > 2)
+	    call error (1, "Image sections only allowed for 1D and 2D images")
+
+        # Decode the section string.
+	call amovki (1, a, 2)
+	call amovki (1, b, 2)
+	call amovki (1, c, 2)
+	do i = 1, IM_NDIM(refim)
+	    b[i] = IM_LEN(refim,i)
+
+        ip = 1
+        while (IS_WHITE(section[ip]))
+            ip = ip + 1
+        if (section[ip] == '[') {
+            ip = ip + 1
+
+	    do i = 1, IM_NDIM(refim) {
+		while (IS_WHITE(section[ip]))
+		    ip = ip + 1
+
+		# Get a:b:c.  Allow notation such as "-*:c"
+		# (or even "-:c") where the step is obviously negative.
+
+		if (ctoi (section, ip, temp) > 0) {                 # a
+		    a[i] = temp
+		    if (section[ip] == ':') {
+			ip = ip + 1
+			if (ctoi (section, ip, b[i]) == 0)             # a:b
+			    goto error_
+		    } else
+			b[i] = a[i]
+		} else if (section[ip] == '-') {                    # -*
+		    temp = a[i]
+		    a[i] = b[i]
+		    b[i] = temp
+		    ip = ip + 1
+		    if (section[ip] == '*')
+			ip = ip + 1
+		} else if (section[ip] == '*')                      # *
+		    ip = ip + 1
+		if (section[ip] == ':') {                           # ..:step
+		    ip = ip + 1
+		    if (ctoi (section, ip, c[i]) == 0)
+			goto error_
+		    else if (c[i] == 0)
+			goto error_
+		}
+		if (a[i] > b[i] && c[i] > 0)
+		    c[i] = -c[i]
+
+		while (IS_WHITE(section[ip]))
+		    ip = ip + 1
+		if (i < IM_NDIM(refim)) {
+		    if (section[ip] != ',')
+			goto error_
+		} else {
+		    if (section[ip] != ']')
+			goto error_
+		}
+		ip = ip + 1
+	    }
+	}
+
+	# In this case make the values be increasing only.
+	do i = 1, IM_NDIM(refim)
+	    if (c[i] < 0) {
+		temp = a[i]
+		a[i] = b[i]
+		b[i] = temp
+		c[i] = -c[i]
+	    }
+
+	# Make the mask.
+	pm = pm_newmask (refim, 16)
+
+	if (flag == INVERT_MASK) {
+	    rop = PIX_SET+PIX_VALUE(1)
+	    call pm_box (pm, 1, 1, IM_LEN(refim,1), IM_LEN(refim,2), rop)
+	    rop = PIX_CLR
+	} else
+	    rop = PIX_SET+PIX_VALUE(1)
+
+	if (c[1] == 1 && c[2] == 1)
+	    call pm_box (pm, a[1], a[2], b[1], b[2], rop)
+
+	else if (c[1] == 1)
+	    for (i=a[2]; i<=b[2]; i=i+c[2])
+		call pm_box (pm, a[1], i, b[1], i, rop)
+
+	else
+	    for (i=a[2]; i<=b[2]; i=i+c[2])
+		for (j=a[1]; j<=b[1]; j=j+c[1])
+		    call pm_point (pm, j, i, rop)
+
+	im = im_pmmapo (pm, refim)
+	dummy = imgl1i (im)		# Force I/O to set header
+	mw = mw_openim (refim)		# Set WCS
+	call mw_saveim (mw, im)
+	call mw_close (mw)
+
+	return (im)
+
+error_
+        call error (1, "Error in image section specification")
+end
+
+
+# XT_PMINVERT -- Invert a pixel mask by changing 0 to 1 and non-zero to zero.
+
+procedure xt_pminvert (pm)
+
+pointer	pm		#I Pixel mask to be inverted
+
+int	i, naxes, axlen[IM_MAXDIM], depth, npix, val
+pointer	sp, v, buf, one
+bool	pm_linenotempty()
+
+begin
+	call pm_gsize (pm, naxes, axlen, depth)
+
+	call smark (sp)
+	call salloc (v, IM_MAXDIM, TY_LONG)
+	call salloc (buf, axlen[1], TY_INT)
+	call salloc (one, 6, TY_INT)
+
+	npix = axlen[1]
+	RLI_LEN(one) = 2
+	RLI_AXLEN(one) = npix
+	Memi[one+3] = 1
+	Memi[one+4] = npix
+	Memi[one+5] = 1
+
+	call amovkl (long(1), Meml[v], IM_MAXDIM)
+	repeat {
+	    if (pm_linenotempty (pm, Meml[v])) {
+		call pmglpi (pm, Meml[v], Memi[buf], 0, npix, 0)
+		do i = 0, npix-1 {
+		    val = Memi[buf+i]
+		    if (val == 0)
+			Memi[buf+i] = 1
+		    else
+			Memi[buf+i] = 0
+		}
+		call pmplpi (pm, Meml[v], Memi[buf], 0, npix, PIX_SRC)
+	    } else
+		call pmplri (pm, Meml[v], Memi[one], 0, npix, PIX_SRC)
+	    
+	    do i = 2, naxes {
+		Meml[v+i-1] = Meml[v+i-1] + 1
+		if (Meml[v+i-1] <= axlen[i])
+		    break
+		else if (i < naxes)
+		    Meml[v+i-1] = 1
+	    }
+	} until (Meml[v+naxes-1] > axlen[naxes])
+
+	call sfree (sp)
+end
+
+
+# XT_MATCH -- Set the pixel mask to match the reference image.
+# This matches sizes and physical coordinates and allows the
+# original mask to be smaller or larger than the reference image.
+# Subsequent use of the pixel mask can then work in the logical
+# coordinates of the reference image.  The mask values are the maximum
+# of the mask values which overlap each reference image pixel.
+# A null input returns a null output.
+
+procedure xt_match (im, refim, match)
+
+pointer	im			#U Pixel mask image pointer
+pointer	refim			#I Reference image pointer
+int	match			#I Match by physical coordinates?
+
+int	i, j, k, l, i1, i2, j1, j2, nc, nl, ncpm, nlpm, nx, val
+double	x1, x2, y1, y2, lt[6], lt1[6], lt2[6]
+long	vold[IM_MAXDIM], vnew[IM_MAXDIM]
+pointer	pm, pmnew, imnew, mw, ctx, cty, bufref, bufpm
+
+int	imstati()
+pointer	pm_open(), mw_openim(), im_pmmapo(), imgl1i(), mw_sctran()
+bool	pm_empty(), pm_linenotempty()
+errchk	pm_open, mw_openim, im_pmmapo
+
+begin
+	if (im == NULL)
+	    return
+
+	# Set sizes.
+	nc = IM_LEN(refim,1)
+	nl = IM_LEN(refim,2)
+	ncpm = IM_LEN(im,1)
+	nlpm = IM_LEN(im,2)
+
+	# If the mask is empty and the sizes are the same then it does not
+	# matter if the two are actually matched in physical coordinates.
+	pm = imstati (im, IM_PMDES)
+	if (pm_empty(pm) && nc == ncpm && nl == nlpm)
+	    return
+
+	# Compute transformation between reference (logical) coordinates
+	# and mask (physical) coordinates if desired.
+
+	mw = mw_openim (im)
+	call mw_gltermd (mw, lt, lt[5], 2)
+	call mw_close (mw)
+
+	if (match == YES) {
+	    mw = mw_openim (refim)
+	    call mw_gltermd (mw, lt2, lt2[5], 2)
+	    call mw_close (mw)
+	} else
+	    call amovd (lt, lt2, 6)
+
+	# Combine lterms.
+	call mw_invertd (lt, lt1, 2)
+	call mw_mmuld (lt1, lt2, lt, 2)
+	call mw_vmuld (lt, lt[5], lt[5], 2)
+	lt[5] = lt2[5] - lt[5]
+	lt[6] = lt2[6] - lt[6]
+	do i = 1, 6
+	    lt[i] = nint (1D6 * (lt[i]-int(lt[i]))) / 1D6 + int(lt[i])
+
+	# Check for a rotation.  For now don't allow any rotation.
+	if (lt[2] != 0. || lt[3] != 0.)
+	    call error (1, "Image and mask have a relative rotation")
+	
+	# Check for an exact match.
+	if (lt[1] == 1D0 && lt[4] == 1D0 && lt[5] == 0D0 && lt[6] == 0D0 &&
+	    nc == ncpm && nl == nlpm)
+	    return
+
+	# Set reference to mask coordinates.
+	mw = mw_openim (im)
+	call mw_sltermd (mw, lt, lt[5], 2)
+	ctx = mw_sctran (mw, "logical", "physical", 1)
+	cty = mw_sctran (mw, "logical", "physical", 2)
+
+	# Create a new pixel mask of the required size and offset.
+	# Do dummy image I/O to set the header.
+	pmnew = pm_open (NULL)
+	call pm_ssize (pmnew, 2, IM_LEN(refim,1), 27)
+	imnew = im_pmmapo (pmnew, NULL)
+	bufref = imgl1i (imnew)
+
+	# Compute region of mask overlapping the reference image.
+	call mw_ctrand (ctx, 1-0.5D0, x1, 1)
+	call mw_ctrand (ctx, nc+0.5D0, x2, 1)
+	i1 = max (1, nint(min(x1,x2)+1D-5))
+	i2 = min (ncpm, nint(max(x1,x2)-1D-5))
+	call mw_ctrand (cty, 1-0.5D0, y1, 1)
+	call mw_ctrand (cty, nl+0.5D0, y2, 1)
+	j1 = max (1, nint(min(y1,y2)+1D-5))
+	j2 = min (nlpm, nint(max(y1,y2)-1D-5))
+
+	# Set the new mask values to the maximum of all mask values falling
+	# within each reference pixel in the overlap region.
+	if (i1 <= i2 && j1 <= j2) {
+	    nx = i2 - i1 + 1
+	    vold[1] = i1
+	    vnew[1] = 1
+
+	    # If the scales are the same then it is just a problem of
+	    # padding.  In this case use range lists for speed.
+	    if (lt[1] == 1D0 && lt[4] == 1D0) {
+		call malloc (bufpm, 3+3*nc, TY_INT)
+		k = nint (lt[5])
+		l = nint (lt[6])
+		do j = max(1-l,j1), min(nl-l,j2) {
+		    vold[2] = j
+		    call pmglri (pm, vold, Memi[bufpm], 0, nc, PIX_SRC)
+		    if (k != 0) {
+			bufref = bufpm
+			do i = 2, Memi[bufpm] {
+			    bufref = bufref + 3
+			    Memi[bufref] = Memi[bufref] + k
+			}
+		    }
+		    vnew[2] = j + l
+		    call pmplri (pmnew, vnew, Memi[bufpm], 0, nc, PIX_SRC)
+		}
+		bufref = NULL
+
+	    # Do all the geometry and pixel size matching.  This can
+	    # be slow.
+	    } else {
+		call malloc (bufpm, nx, TY_INT)
+		call malloc (bufref, nc, TY_INT)
+		do j = 1, nl {
+		    call mw_ctrand (cty, j-0.5D0, y1, 1)
+		    call mw_ctrand (cty, j+0.5D0, y2, 1)
+		    j1 = max (1, nint(min(y1,y2)+1D-5))
+		    j2 = min (nlpm, nint(max(y1,y2)-1D-5))
+		    if (j2 < j1)
+			next
+
+		    vnew[2] = j
+		    call aclri (Memi[bufref], nc)
+		    do l = j1, j2 {
+			vold[2] = l
+			if (!pm_linenotempty (pm, vold))
+			    next
+			call pmglpi (pm, vold, Memi[bufpm], 0, nx, 0)
+			do i = 1, nc {
+			    call mw_ctrand (ctx, i-0.5D0, x1, 1)
+			    call mw_ctrand (ctx, i+0.5D0, x2, 1)
+			    i1 = max (1, nint(min(x1,x2)+1D-5))
+			    i2 = min (ncpm, nint(max(x1,x2)-1D-5))
+			    if (i2 < i1)
+				next
+			    val = Memi[bufref+i-1]
+			    do k = i1-vold[1], i2-vold[1]
+				val = max (val, Memi[bufpm+k])
+			    Memi[bufref+i-1] = val
+			}
+		    }
+		    call pmplpi (pmnew, vnew, Memi[bufref], 0, nc, PIX_SRC)
+		}
+	    }
+	    call mfree (bufref, TY_INT)
+	    call mfree (bufpm, TY_INT)
+	}
+
+	call mw_close (mw)
+	call xt_pmunmap (im)
+	im = imnew
+	call imseti (im, IM_PMDES, pmnew)
+end
diff --git a/pkg/xtools/fixpix/ytfixpix.x b/pkg/xtools/fixpix/ytfixpix.x
new file mode 100644
index 00000000..e93b4c07
--- /dev/null
+++ b/pkg/xtools/fixpix/ytfixpix.x
@@ -0,0 +1,281 @@
+include	<imhdr.h>
+include	<imset.h>
+include	<pmset.h>
+include	"xtfixpix.h"
+
+# This version uses an internal copy of the input mask rather than modifying
+# the input mask.
+
+
+# XT_FPINIT -- Initialize FIXPIX data structure.
+# If the mask is null or empty a null pointer is returned.
+# If the mask is not empty the mask is examined for bad pixels requiring
+# column interpolation.  The columns and interpolation endpoints are
+# recorded.  Note that line interpolation does not need to be mapped since
+# this can be done efficiently as the reference image is accessed line by
+# line.
+
+pointer procedure yt_fpinit (pmin, lvalin, cvalin)
+
+pointer	pmin			#I Pixel mask
+int	lvalin			#I Input line interpolation code
+int	cvalin			#I Input column interpolation code
+
+int	i, j, k, l, n, nc, nl, l1, l2, lmin, lmax, ncols, lval, cval, ncompress
+short	val
+long	v[IM_MAXDIM]
+pointer	pm, fp, ptr, col, pl1, pl2
+pointer	sp, buf, cols
+
+bool	pm_empty()
+pointer	pm_newcopy()
+errchk	pmglrs, pmplrs
+
+begin
+	# Check for empty mask.
+	if (pmin == NULL)
+	    return (NULL)
+	if (pm_empty (pmin))
+	    return (NULL)
+
+	# Make an internal copy of the mask.
+	pm = pm_newcopy (pmin)
+
+	# Get mask size.
+	call pm_gsize (pm, i, v, j)
+	nc = v[1]
+	nl = v[2]
+
+	# Allocate memory and data structure.
+	call smark (sp)
+	call salloc (buf, 3*max(nc, nl), TY_SHORT) 
+	call salloc (cols, nc, TY_SHORT)
+	call calloc (fp, FP_LEN, TY_STRUCT)
+
+	# Set the mask codes.  Go through the mask and change any mask codes
+	# that match the input mask code to the output mask code (if they are
+	# different).  This is done to move the mask codes to a range that
+	# won't conflict with the length values.  For any other code replace
+	# the value by the length of the bad region along the line.  This
+	# value will be used in comparison to the length along the column for
+	# setting the interpolation for the narrower dimension.
+
+	if ((IS_INDEFI(lvalin)||lvalin<1) && (IS_INDEFI(cvalin)||cvalin<1)) {
+	    lval = FP_LDEF
+	    cval = FP_CDEF
+	} else if (IS_INDEFI(lvalin) || lvalin < 1) {
+	    lval = FP_LDEF
+	    cval = mod (cvalin - 1, nc) + 1
+	    if (lval == cval)
+		lval = FP_CDEF
+	} else if (IS_INDEFI(cvalin) || cvalin < 1) {
+	    lval = mod (lvalin - 1, nc) + 1
+	    cval = FP_CDEF
+	    if (cval == lval)
+		cval = FP_LDEF
+	} else if (lvalin != cvalin) {
+	    lval = mod (lvalin - 1, nc) + 1
+	    cval = mod (cvalin - 1, nc) + 1
+	} else {
+	    call mfree (fp, TY_STRUCT)
+	    call sfree (sp)
+	    call error (1, "Interpolation codes cannot be the same")
+	}
+	call yt_fpsinterp (pmin, pm, nc, nl, v, Mems[buf], lvalin, cvalin,
+	    lval, cval)
+
+	# Go through and check if there is any need for column interpolation;
+	# i.e. are there any mask values different from the line interpolation.
+
+	call aclrs (Mems[cols], nc)
+	call amovkl (long(1), v, IM_MAXDIM)
+	do l = 1, nl {
+	    v[2] = l
+	    call pmglrs (pm, v, Mems[buf], 0, nc, 0)
+	    ptr = buf + 3
+	    do i = 2, Mems[buf] {
+		val = Mems[ptr+2]
+		if (val != lval) {
+		    val = 1
+		    n = Mems[ptr+1] 
+		    call amovks (val, Mems[cols+Mems[ptr]-1], n)
+		}
+		ptr = ptr + 3
+	    }
+	}
+	n = 0
+	do i = 1, nc
+	    if (Mems[cols+i-1] != 0)
+		n = n + 1
+
+	# If there are mask codes for either column interpolation or
+	# interpolation lengths along lines to compare against column
+	# interpolation check the interpolation length against the
+	# column and set the line interpolation endpoints to use.
+	# compute the minimum and maximum lines that are endpoints
+	# to restrict the random access pass that will be needed to
+	# get the endpoint values.
+
+	if (n > 0) {
+	    n = n + 10
+	    call malloc (col, n, TY_INT)
+	    call malloc (pl1, n, TY_INT)
+	    call malloc (pl2, n, TY_INT)
+	    ncols = 0
+	    lmin = nl
+	    lmax = 0
+	    ncompress = 0
+	    do i = 1, nc {
+		if (Mems[cols+i-1] == 0)
+		    next
+		v[1] = i
+		do l = 1, nl {
+		    v[2] = l
+		    call pmglps (pm, v, Mems[buf+l-1], 0, 1, 0)
+		}
+		for (l1=1; l1<=nl && Mems[buf+l1-1]==0; l1=l1+1)
+		    ;
+		while (l1 <= nl) {
+		    l1 = l1 - 1
+		    for (l2=l1+1; l2<=nl && Mems[buf+l2-1]!=0; l2=l2+1)
+			;
+		    j = 0
+		    k = nc + l2 - l1 - 1
+		    do l = l1+1, l2-1 {
+			val = Mems[buf+l-1]
+			if (val == cval)
+			    j = j + 1
+			else if (val > nc) {
+			    if (val > k) {
+				j = j + 1
+				val = cval
+			    } else
+				val = lval
+			    v[2] = l
+			    call pmplps (pm, v, val, 0, 1, PIX_SRC)
+			    ncompress = ncompress + 1
+			}
+		    }
+		    if (ncompress > 100) {
+			call pm_compress (pm)
+			ncompress = 0
+		    }
+		    if (j > 0) {
+			if (ncols == n) {
+			    n = n + 10
+			    call realloc (col, n, TY_INT)
+			    call realloc (pl1, n, TY_INT)
+			    call realloc (pl2, n, TY_INT)
+			}
+			j = 1 + l1 - 1
+			k = 1 + l2 - 1
+			lmin = min (lmin, j, k)
+			lmax = max (lmax, j, k)
+			Memi[col+ncols] = i
+			Memi[pl1+ncols] = j
+			Memi[pl2+ncols] = k
+			ncols = ncols + 1
+		    }
+		    for (l1=l2+1; l1<=nl && Mems[buf+l1-1]==0; l1=l1+1)
+			;
+		}
+	    }
+
+	    FP_LMIN(fp) = lmin
+	    FP_LMAX(fp) = lmax
+	    FP_NCOLS(fp) = ncols
+	    FP_PCOL(fp) = col
+	    FP_PL1(fp) = pl1
+	    FP_PL2(fp) = pl2
+	}
+
+	FP_PM(fp) = pm
+	FP_LVAL(fp) = lval
+	FP_CVAL(fp) = cval
+
+	call sfree (sp)
+	return (fp)
+end
+
+
+# XT_SINTERP -- Set length of line interpolation regions.
+# The mask values are set to the length of any column interpolation
+# plus an offset leaving any line and column interpolation codes
+# unchanged.  These values will be used in a second pass to compare
+# to the lengths of line interpolation and then the mask values will
+# be reset to one of the line or column interpolation codes based on
+# the minimum distance.
+
+procedure yt_fpsinterp (pmin, pm, nc, nl, v, data, lvalin, cvalin,
+	lvalout, cvalout)
+
+pointer	pmin		#I Input pixel mask
+pointer	pm		#I Modified pixel mask
+int	nc, nl		#I Mask size
+long	v[ARB]		#I Coordinate vector
+short	data[ARB]	#I Data buffer
+int	lvalin		#I Input line interpolation code
+int	cvalin		#I Input column interpolation code
+int	lvalout		#I Output line interpolation code
+int	cvalout		#I Output column interpolation code
+
+int	c, l, c1, c2, val
+bool	pm_linenotempty()
+
+begin
+	call amovkl (long(1), v, IM_MAXDIM)
+	do l = 1, nl {
+	    v[2] = l
+	    if (!pm_linenotempty (pmin, v))
+		next
+	    
+	    call pmglps (pmin, v, data, 0, nc, 0)
+
+	    for (c1=1; c1<=nc && data[c1]==0; c1=c1+1)
+		;
+	    while (c1 <= nc) {
+		for (c2=c1+1; c2<=nc && data[c2]!=0; c2=c2+1)
+		    ;
+		c2 = c2 - 1
+		do c = c1, c2 {
+		    val = data[c]
+		    if (val == lvalin) {
+			if (lvalin != lvalout)
+			    data[c] = lvalout
+		    } else if (val == cvalin) {
+			if (cvalin != cvalout)
+			    data[c] = cvalout
+		    } else {
+			data[c] = nc + c2 - c1 + 1
+		    }
+		}
+		for (c1=c2+2; c1<=nc && data[c1]==0; c1=c1+1)
+		    ;
+	    }
+
+	    call pmplps (pm, v, data, 0, nc, PIX_SRC)
+	}
+end
+
+
+# XT_FPFREE -- Free FIXPIX data structures.
+
+procedure yt_fpfree (fp)
+
+pointer	fp		#U FIXPIX data structure
+
+begin
+	if (fp == NULL)
+	    return
+	call mfree (FP_PCOL(fp), TY_INT)
+	call mfree (FP_PL1(fp), TY_INT)
+	call mfree (FP_PL2(fp), TY_INT)
+	if (FP_PV1(fp) != NULL)
+	    call mfree (FP_PV1(fp), FP_PIXTYPE(fp))
+	if (FP_PV2(fp) != NULL)
+	    call mfree (FP_PV2(fp), FP_PIXTYPE(fp))
+	if (FP_DATA(fp) != NULL)
+	    call mfree (FP_DATA(fp), FP_PIXTYPE(fp))
+	call pm_close (FP_PM(fp))
+	call mfree (fp, TY_STRUCT)
+end
diff --git a/pkg/xtools/fixpix/ytpmmap.x b/pkg/xtools/fixpix/ytpmmap.x
new file mode 100644
index 00000000..e41fb4f8
--- /dev/null
+++ b/pkg/xtools/fixpix/ytpmmap.x
@@ -0,0 +1,961 @@
+include	<mach.h>
+include	<ctype.h>
+include	<error.h>
+include	<imhdr.h>
+include	<imset.h>
+include	<pmset.h>
+include	<mwset.h>
+include	<syserr.h>
+include	<math/iminterp.h>
+
+# Pixel mask matching options.
+define	PM_MATCH	"|logical|physical|world|offset|"
+define	PM_LOGICAL	1		# Match in logical coordinates
+define	PM_PHYSICAL	2		# Match in physical coordinates
+define	PM_WORLD	3		# Match in world coordinates
+define	PM_OFFSET	4		# Match in physical with WCS offset
+
+
+# XT_PMMAP/XT_MAPPM -- Open a pixel mask READ_ONLY.
+#
+# This routine maps multiple types of mask files and designations.
+# It may match the mask coordinates to the reference image based on the
+# physical coordinate system so the mask may be of a different size.
+# The mask name is returned so that the task has the name pointed to by "BPM".
+# A null filename is allowed and returns NULL.
+#
+# Modified to use xt_maskname with the reference image extension name.
+# Minor bug fixes in xt_match.
+
+pointer procedure yt_pmmap (pmname, refim, mname, sz_mname)
+
+char	pmname[ARB]		#I Pixel mask name
+pointer	refim			#I Reference image pointer
+char	mname[ARB]		#O Expanded mask name
+int	sz_mname		#O Size of expanded mask name
+
+pointer	yt_mappm()
+errchk	yt_mappm
+
+begin
+	return (yt_mappm (pmname, refim, "physical", mname, sz_mname))
+end
+
+pointer procedure yt_mappm (pmname, refim, match, mname, sz_mname)
+
+char	pmname[ARB]		#I Pixel mask name
+pointer	refim			#I Reference image pointer
+char	match[ARB]		#I Match by physical coordinates?
+char	mname[ARB]		#O Expanded mask name
+int	sz_mname		#O Size of expanded mask name
+
+int	i, j, flag, nowhite()
+pointer	sp, fname, extname, im, ref, yt_pmmap1()
+bool	streq()
+errchk	yt_pmmap1
+
+begin
+	call smark (sp)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+	call salloc (extname, SZ_FNAME, TY_CHAR)
+
+	im = NULL
+	i = nowhite (pmname, Memc[fname], SZ_FNAME)
+
+	# Process invert flags.  These occur more than once.
+	j = 0; flag = 0
+	for (i=0; Memc[fname+i]!=EOS; i=i+1) {
+	    if (Memc[fname+i] == '^')
+	        flag = flag + 1
+	    else {
+	        Memc[fname+j] = Memc[fname+i]
+		j = j + 1
+	    }
+	}
+	Memc[fname+j] = EOS
+	if (mod (flag, 2) == 0)
+	    flag = 0
+	else
+	    flag = INVERT_MASK
+
+
+	# Resolve keyword references.
+	if (Memc[fname] == '!') {
+	    iferr (call imgstr (refim, Memc[fname+1], Memc[fname], SZ_FNAME))
+		Memc[fname] = EOS
+	} else if (streq (Memc[fname], "BPM")) {
+	    iferr (call imgstr (refim, "BPM", Memc[fname], SZ_FNAME))
+		Memc[fname] = EOS
+	}
+
+	# Resolve other special names.
+	if (streq (Memc[fname], "EMPTY"))
+	    ref = refim
+	else
+	    ref = NULL
+
+	# Create the mask.
+	if (Memc[fname] != EOS) {
+	    iferr (im = yt_pmmap1 (Memc[fname], ref, refim, flag, match)) {
+	        ifnoerr (call imgstr (refim, "extname", Memc[extname],
+		    SZ_FNAME)) {
+		    call xt_maskname (Memc[fname], Memc[extname], READ_ONLY,
+		        Memc[fname], SZ_FNAME)
+		    im = yt_pmmap1 (Memc[fname], ref, refim, flag, match)
+		} else
+		    im = yt_pmmap1 (Memc[fname], ref, refim, flag, match)
+	    }
+	}
+	call strcpy (Memc[fname], mname, sz_mname)
+
+	call sfree (sp)
+	return (im)
+end
+
+
+# XT_PMUNMAP -- Unmap a mask image.
+# Note that the imio pointer may be purely an internal pointer opened
+# with im_pmmapo so we need to free the pl pointer explicitly.
+
+procedure yt_pmunmap (im)
+
+pointer	im			#I IMIO pointer for mask
+
+pointer	pm
+int	imstati()
+
+begin
+	pm = imstati (im, IM_PMDES)
+	call pm_close (pm)
+	call imseti (im, IM_PMDES, NULL)
+	call imunmap (im)
+end
+
+
+# XT_PMMAP1 -- Open a pixel mask READ_ONLY.  The input mask may be
+# a pixel list image, a non-pixel list image, or a text file.
+# Return error if the pixel mask cannot be opened.  For pixel masks
+# or image masks possibly match the WCS.
+
+pointer procedure yt_pmmap1 (pmname, ref, refim, flag, match)
+
+char	pmname[ARB]		#I Pixel mask name
+pointer	ref			#I Reference image for pixel mask
+pointer	refim			#I Reference image for image or text
+int	flag			#I Mask flag
+char	match[ARB]		#I Match by physical coordinates?
+
+int	imstati(),  errcode()
+pointer	im, pm
+pointer	im_pmmap(), yt_pmimmap(), yt_pmtext(), yt_pmsection()
+bool	streq()
+errchk	yt_match
+
+begin
+	im = NULL
+
+	if (streq (pmname, "STDIN"))
+	    im = yt_pmtext (pmname, refim, flag)
+
+	else if (pmname[1] == '[')
+	    im = yt_pmsection (pmname, refim, flag)
+
+	else {
+	    ifnoerr (im = im_pmmap (pmname, READ_ONLY, ref)) {
+		call yt_match (im, refim, match)
+		if (flag == INVERT_MASK) {
+		    pm = imstati (im, IM_PMDES)
+		    call yt_pminvert (pm)
+		    call imseti (im, IM_PMDES, pm)
+		}
+	    } else {
+		switch (errcode()) {
+		case SYS_IKIOPEN, SYS_FOPNNEXFIL, SYS_PLBADSAVEF, SYS_FOPEN:
+		    ifnoerr (im = yt_pmimmap (pmname, refim, flag))
+			call yt_match (im, refim, match)
+		    else {
+			switch (errcode()) {
+			case SYS_IKIOPEN:
+			    im = yt_pmtext (pmname, refim, flag)
+			default:
+			    call erract (EA_ERROR)
+			}
+		    }
+		default:
+		    call erract (EA_ERROR)
+		}
+	    }
+	}
+
+	return (im)
+end
+
+
+# XT_PMIMMAP -- Open a pixel mask from a non-pixel list image.
+# Return error if the image cannot be opened.
+
+pointer procedure yt_pmimmap (pmname, refim, flag)
+
+char	pmname[ARB]		#I Image name
+pointer	refim			#I Reference image pointer
+int	flag			#I Mask flag
+
+int	i, ndim, npix, rop, val
+pointer	sp, v1, v2, im_in, im_out, pm, mw, data
+
+int	imstati(), imgnli()
+pointer immap(), pm_newmask(), im_pmmapo(), imgl1i(), mw_openim()
+errchk	immap, mw_openim, im_pmmapo
+
+begin
+	call smark (sp)
+	call salloc (v1, IM_MAXDIM, TY_LONG)
+	call salloc (v2, IM_MAXDIM, TY_LONG)
+
+	call amovkl (long(1), Meml[v1], IM_MAXDIM)
+	call amovkl (long(1), Meml[v2], IM_MAXDIM)
+
+	im_in = immap (pmname, READ_ONLY, 0)
+	pm = imstati (im_in, IM_PMDES)
+	if (pm != NULL)
+	    return (im_in)
+	pm = pm_newmask (im_in, 16)
+
+	ndim = IM_NDIM(im_in)
+	npix = IM_LEN(im_in,1)
+
+	if (flag == INVERT_MASK)
+	    rop = PIX_NOT(PIX_SRC)
+	else
+	    rop = PIX_SRC
+
+	while (imgnli (im_in, data, Meml[v1]) != EOF) {
+	    if (flag == INVERT_MASK) {
+		do i = 0, npix-1 {
+		    val = Memi[data+i]
+		    if (val <= 0)
+		       Memi[data+i] = 1
+		    else
+		       Memi[data+i] = 0
+		}
+	    } else {
+		do i = 0, npix-1 {
+		    val = Memi[data+i]
+		    if (val < 0)
+		       Memi[data+i] = 0
+		}
+	    }
+	    call pmplpi (pm, Meml[v2], Memi[data], 0, npix, rop)
+	    call amovl (Meml[v1], Meml[v2], ndim)
+	}
+
+	im_out = im_pmmapo (pm, im_in)
+	data = imgl1i (im_out)		# Force I/O to set header
+	mw = mw_openim (im_in)		# Set WCS
+	call mw_saveim (mw, im_out)
+	call mw_close (mw)
+
+	#call imunmap (im_in)
+	call yt_pmunmap (im_in)
+	call sfree (sp)
+	return (im_out)
+end
+
+
+# XT_PMTEXT -- Create a pixel mask from a text file of rectangles.
+# Return error if the file cannot be opened.
+# This routine only applies to the first 2D plane.
+
+pointer procedure yt_pmtext (pmname, refim, flag)
+
+char	pmname[ARB]		#I Image name
+pointer	refim			#I Reference image pointer
+int	flag			#I Mask flag
+
+int	fd, nc, nl, c1, c2, l1, l2, nc1, nl1, rop
+pointer	pm, im, mw, dummy
+
+int	open(), fscan(), nscan()
+pointer	pm_newmask(), im_pmmapo(), imgl1i(), mw_openim()
+errchk	open,im_pmmapo
+
+begin
+	fd = open (pmname, READ_ONLY, TEXT_FILE)
+	pm = pm_newmask (refim, 16)
+
+	nc = IM_LEN(refim,1)
+	nl = IM_LEN(refim,2)
+
+	if (flag == INVERT_MASK)
+	    call pl_box (pm, 1, 1, nc, nl, PIX_SET+PIX_VALUE(1))
+
+	while (fscan (fd) != EOF) {
+	    call gargi (c1)
+	    call gargi (c2)
+	    call gargi (l1)
+	    call gargi (l2)
+	    if (nscan() != 4) {
+		if (nscan() == 2) {
+		    l1 = c2
+		    c2 = c1
+		    l2 = l1
+		} else
+		    next
+	    }
+
+	    c1 = max (1, c1)
+	    c2 = min (nc, c2)
+	    l1 = max (1, l1)
+	    l2 = min (nl, l2)
+	    nc1 = c2 - c1 + 1
+	    nl1 = l2 - l1 + 1
+	    if (nc1 < 1 || nl1 < 1)
+		next
+
+	    # Select mask value based on shape of rectangle.
+	    if (flag == INVERT_MASK)
+		rop = PIX_CLR
+	    else if (nc1 <= nl1)
+		rop = PIX_SET+PIX_VALUE(2)
+	    else
+		rop = PIX_SET+PIX_VALUE(3)
+
+	    # Set mask rectangle.
+	    call pm_box (pm, c1, l1, c2, l2, rop)
+	}
+
+	call close (fd)
+	im = im_pmmapo (pm, refim)
+	dummy = imgl1i (im)		# Force I/O to set header
+	mw = mw_openim (refim)		# Set WCS
+	call mw_saveim (mw, im)
+	call mw_close (mw)
+
+	return (im)
+end
+
+
+# XT_PMSECTION -- Create a pixel mask from an image section.
+# This only applies the mask to the first plane of the image.
+
+pointer procedure yt_pmsection (section, refim, flag)
+
+char	section[ARB]		#I Image section
+pointer	refim			#I Reference image pointer
+int	flag			#I Mask flag
+
+int	i, j, ip, temp, a[2], b[2], c[2], rop, ctoi()
+pointer	pm, im, mw, dummy, pm_newmask(), im_pmmapo(), imgl1i(), mw_openim()
+errchk	im_pmmapo
+define  error_  99
+
+begin
+	# This is currently only for 1D and 2D images.
+	if (IM_NDIM(refim) > 2)
+	    call error (1, "Image sections only allowed for 1D and 2D images")
+
+        # Decode the section string.
+	call amovki (1, a, 2)
+	call amovki (1, b, 2)
+	call amovki (1, c, 2)
+	do i = 1, IM_NDIM(refim)
+	    b[i] = IM_LEN(refim,i)
+
+        ip = 1
+        while (IS_WHITE(section[ip]))
+            ip = ip + 1
+        if (section[ip] == '[') {
+            ip = ip + 1
+
+	    do i = 1, IM_NDIM(refim) {
+		while (IS_WHITE(section[ip]))
+		    ip = ip + 1
+
+		# Get a:b:c.  Allow notation such as "-*:c"
+		# (or even "-:c") where the step is obviously negative.
+
+		if (ctoi (section, ip, temp) > 0) {                 # a
+		    a[i] = temp
+		    if (section[ip] == ':') {
+			ip = ip + 1
+			if (ctoi (section, ip, b[i]) == 0)             # a:b
+			    goto error_
+		    } else
+			b[i] = a[i]
+		} else if (section[ip] == '-') {                    # -*
+		    temp = a[i]
+		    a[i] = b[i]
+		    b[i] = temp
+		    ip = ip + 1
+		    if (section[ip] == '*')
+			ip = ip + 1
+		} else if (section[ip] == '*')                      # *
+		    ip = ip + 1
+		if (section[ip] == ':') {                           # ..:step
+		    ip = ip + 1
+		    if (ctoi (section, ip, c[i]) == 0)
+			goto error_
+		    else if (c[i] == 0)
+			goto error_
+		}
+		if (a[i] > b[i] && c[i] > 0)
+		    c[i] = -c[i]
+
+		while (IS_WHITE(section[ip]))
+		    ip = ip + 1
+		if (i < IM_NDIM(refim)) {
+		    if (section[ip] != ',')
+			goto error_
+		} else {
+		    if (section[ip] != ']')
+			goto error_
+		}
+		ip = ip + 1
+	    }
+	}
+
+	# In this case make the values be increasing only.
+	do i = 1, IM_NDIM(refim)
+	    if (c[i] < 0) {
+		temp = a[i]
+		a[i] = b[i]
+		b[i] = temp
+		c[i] = -c[i]
+	    }
+
+	# Make the mask.
+	pm = pm_newmask (refim, 16)
+
+	if (flag == INVERT_MASK) {
+	    rop = PIX_SET+PIX_VALUE(1)
+	    call pm_box (pm, 1, 1, IM_LEN(refim,1), IM_LEN(refim,2), rop)
+	    rop = PIX_CLR
+	} else
+	    rop = PIX_SET+PIX_VALUE(1)
+
+	if (c[1] == 1 && c[2] == 1)
+	    call pm_box (pm, a[1], a[2], b[1], b[2], rop)
+
+	else if (c[1] == 1)
+	    for (i=a[2]; i<=b[2]; i=i+c[2])
+		call pm_box (pm, a[1], i, b[1], i, rop)
+
+	else
+	    for (i=a[2]; i<=b[2]; i=i+c[2])
+		for (j=a[1]; j<=b[1]; j=j+c[1])
+		    call pm_point (pm, j, i, rop)
+
+	im = im_pmmapo (pm, refim)
+	dummy = imgl1i (im)		# Force I/O to set header
+	mw = mw_openim (refim)		# Set WCS
+	call mw_saveim (mw, im)
+	call mw_close (mw)
+
+	return (im)
+
+error_
+        call error (1, "Error in image section specification")
+end
+
+
+# XT_PMINVERT -- Invert a pixel mask by changing 0 to 1 and non-zero to zero.
+
+procedure yt_pminvert (pm)
+
+pointer	pm		#I Pixel mask to be inverted
+
+int	i, naxes, axlen[IM_MAXDIM], depth, npix, val
+pointer	sp, v, buf, one
+bool	pm_linenotempty()
+
+begin
+	call pm_gsize (pm, naxes, axlen, depth)
+
+	call smark (sp)
+	call salloc (v, IM_MAXDIM, TY_LONG)
+	call salloc (buf, axlen[1], TY_INT)
+	call salloc (one, 6, TY_INT)
+
+	npix = axlen[1]
+	RLI_LEN(one) = 2
+	RLI_AXLEN(one) = npix
+	Memi[one+3] = 1
+	Memi[one+4] = npix
+	Memi[one+5] = 1
+
+	call amovkl (long(1), Meml[v], IM_MAXDIM)
+	repeat {
+	    if (pm_linenotempty (pm, Meml[v])) {
+		call pmglpi (pm, Meml[v], Memi[buf], 0, npix, 0)
+		do i = 0, npix-1 {
+		    val = Memi[buf+i]
+		    if (val == 0)
+			Memi[buf+i] = 1
+		    else
+			Memi[buf+i] = 0
+		}
+		call pmplpi (pm, Meml[v], Memi[buf], 0, npix, PIX_SRC)
+	    } else
+		call pmplri (pm, Meml[v], Memi[one], 0, npix, PIX_SRC)
+	    
+	    do i = 2, naxes {
+		Meml[v+i-1] = Meml[v+i-1] + 1
+		if (Meml[v+i-1] <= axlen[i])
+		    break
+		else if (i < naxes)
+		    Meml[v+i-1] = 1
+	    }
+	} until (Meml[v+naxes-1] > axlen[naxes])
+
+	call sfree (sp)
+end
+
+
+# XT_MATCH -- Set the pixel mask to match the reference image.
+# This matches sizes and possibly the physical coordinates and allows the
+# original mask to be smaller or larger than the reference image.
+# Subsequent use of the pixel mask can then work in the logical
+# coordinates of the reference image.  The mask values are the maximum
+# of the mask values which overlap each reference image pixel.
+# A null input returns a null output.
+
+procedure yt_match (im, refim, match)
+
+pointer	im			#U Pixel mask image pointer
+pointer	refim			#I Reference image pointer
+char	match[ARB]		#I Match by physical coordinates?
+
+int	i, j, k, l, i1, i2, j1, j2, nc, nl, ncpm, nlpm, nx, val
+int	pmmatch, maxmaskval
+double	x1, x2, y1, y2, lt[6], lt1[6], lt2[6]
+long	vold[IM_MAXDIM], vnew[IM_MAXDIM]
+pointer	str, pm, pmnew, imnew, mw, ctx, cty, bufref, bufpm
+
+int	imstati(), strdic(), envfind(), nscan()
+pointer	pm_open(), mw_openim(), im_pmmapo(), imgl1i(), mw_sctran()
+bool	pm_empty(), pm_linenotempty()
+errchk	yt_match_world, pm_open, mw_openim, im_pmmapo
+
+begin
+	if (im == NULL)
+	    return
+
+	# Set sizes.
+	nc = IM_LEN(refim,1)
+	nl = IM_LEN(refim,2)
+	ncpm = IM_LEN(im,1)
+	nlpm = IM_LEN(im,2)
+
+	# If the mask is empty and the sizes are the same then it does not
+	# matter if the two are actually matched in physical coordinates.
+	pm = imstati (im, IM_PMDES)
+	if (pm_empty(pm) && nc == ncpm && nl == nlpm)
+	    return
+
+	# Set match type.
+	call malloc (str, SZ_FNAME, TY_CHAR)
+	call sscan (match)
+	call gargwrd (Memc[str], SZ_FNAME); call gargi (maxmaskval)
+	if (nscan() == 1)
+	    maxmaskval = 1
+	pmmatch = strdic (Memc[str], Memc[str], SZ_FNAME, PM_MATCH)
+	if (pmmatch == 0 && match[1] != EOS) {
+	    if (envfind (match, Memc[str], SZ_FNAME) > 0) {
+		call sscan (Memc[str])
+		call gargwrd (Memc[str], SZ_FNAME); call gargi (maxmaskval)
+		if (nscan() == 1)
+		    maxmaskval = 1
+		pmmatch = strdic (Memc[str], Memc[str], SZ_FNAME, PM_MATCH)
+	    } else
+	        pmmatch = PM_LOGICAL
+	} else {
+	    if (envfind ("pmmatch", Memc[str], SZ_FNAME) > 0) {
+		call sscan (Memc[str])
+		call gargwrd (Memc[str], SZ_FNAME); call gargi (maxmaskval)
+		if (nscan() == 1)
+		    maxmaskval = 1
+		pmmatch = strdic (Memc[str], Memc[str], SZ_FNAME, PM_MATCH)
+	    }
+	}
+	call mfree (str, TY_CHAR)
+	if (pmmatch == 0)
+	    call error (1, "Unknown or invalid pixel mask matching option")
+
+	if (pmmatch == PM_WORLD) {
+	    call yt_match_world (im, refim, maxmaskval)
+	    return
+	}
+
+	# Compute transformation between reference (logical) coordinates
+	# and mask (physical) coordinates.  Apply a world coordinate
+	# offset if desired.
+
+	mw = mw_openim (im)
+	if (pmmatch == PM_OFFSET) {
+	    call mw_gwtermd (mw, lt[5], lt1, lt, 2)
+	    ctx = mw_sctran (mw, "world", "physical", 0)
+	    call mw_ctrand (ctx, lt1, lt1[5], 2)
+	} else
+	    call aclrd (lt1[5], 2)
+	call mw_gltermd (mw, lt, lt[5], 2)
+	call mw_close (mw)
+
+	if (pmmatch == PM_LOGICAL)
+	    call amovd (lt, lt2, 6)
+	else {
+	    mw = mw_openim (refim)
+	    if (pmmatch == PM_OFFSET) {
+		ctx = mw_sctran (mw, "world", "physical", 0)
+		call mw_ctrand (ctx, lt1, lt1[3], 2)
+		lt1[5] = nint (lt1[5] - lt1[3])
+		lt1[6] = nint (lt1[6] - lt1[4])
+	    }
+	    call mw_gltermd (mw, lt2, lt2[5], 2)
+	    lt2[5] = lt2[5] - lt1[5]
+	    lt2[6] = lt2[6] - lt1[6]
+	    call mw_close (mw)
+	}
+
+	# Combine lterms.
+	call mw_invertd (lt, lt1, 2)
+	call mw_mmuld (lt1, lt2, lt, 2)
+	call mw_vmuld (lt, lt[5], lt[5], 2)
+	lt[5] = lt2[5] - lt[5]
+	lt[6] = lt2[6] - lt[6]
+	do i = 1, 6
+	    lt[i] = nint (1D6 * (lt[i]-int(lt[i]))) / 1D6 + int(lt[i])
+
+	# Check for a rotation.  For now don't allow any rotation.
+	if (lt[2] != 0. || lt[3] != 0.)
+	    call error (1, "Image and mask have a relative rotation")
+	
+	# Check for an exact match.
+	if (lt[1] == 1D0 && lt[4] == 1D0 && lt[5] == 0D0 && lt[6] == 0D0 &&
+	    nc == ncpm && nl == nlpm)
+	    return
+
+	# Set reference to mask coordinates.
+	mw = mw_openim (im)
+	call mw_sltermd (mw, lt, lt[5], 2)
+	ctx = mw_sctran (mw, "logical", "physical", 1)
+	cty = mw_sctran (mw, "logical", "physical", 2)
+
+	# Create a new pixel mask of the required size and offset.
+	# Do dummy image I/O to set the header.
+	pmnew = pm_open (NULL)
+	call pm_ssize (pmnew, 2, IM_LEN(refim,1), 27)
+	imnew = im_pmmapo (pmnew, NULL)
+	bufref = imgl1i (imnew)
+
+	# Compute region of mask overlapping the reference image.
+	call mw_ctrand (ctx, 1-0.5D0, x1, 1)
+	call mw_ctrand (ctx, nc+0.5D0, x2, 1)
+	i1 = max (1, nint(min(x1,x2)+1D-5))
+	i2 = min (ncpm, nint(max(x1,x2)-1D-5))
+	call mw_ctrand (cty, 1-0.5D0, y1, 1)
+	call mw_ctrand (cty, nl+0.5D0, y2, 1)
+	j1 = max (1, nint(min(y1,y2)+1D-5))
+	j2 = min (nlpm, nint(max(y1,y2)-1D-5))
+
+	# Set the new mask values to the maximum of all mask values falling
+	# within each reference pixel in the overlap region.
+	if (i1 <= i2 && j1 <= j2) {
+	    nx = i2 - i1 + 1
+	    vold[1] = i1
+	    vnew[1] = 1
+
+	    # If the scales are the same then it is just a problem of
+	    # padding.  In this case use range lists for speed.
+	    if (lt[1] == 1D0 && lt[4] == 1D0) {
+		call malloc (bufpm, 3+3*nc, TY_INT)
+		k = nint (lt[5])
+		l = nint (lt[6])
+		do j = max(1-l,j1), min(nl-l,j2) {
+		    vold[2] = j
+		    call plglri (pm, vold, Memi[bufpm], 0, nc, PIX_SRC)
+		    if (k != 0) {
+			bufref = bufpm
+			do i = 2, Memi[bufpm] {
+			    bufref = bufref + 3
+			    Memi[bufref] = Memi[bufref] + k
+			}
+		    }
+		    vnew[2] = j + l
+		    call pmplri (pmnew, vnew, Memi[bufpm], 0, nc, PIX_SRC)
+		}
+		bufref = NULL
+
+	    # Do all the geometry and pixel size matching.  This can
+	    # be slow.
+	    } else {
+		call malloc (bufpm, nx, TY_INT)
+		call malloc (bufref, nc, TY_INT)
+		do j = 1, nl {
+		    call mw_ctrand (cty, j-0.5D0, y1, 1)
+		    call mw_ctrand (cty, j+0.5D0, y2, 1)
+		    j1 = max (1, nint(min(y1,y2)+1D-5))
+		    j2 = min (nlpm, nint(max(y1,y2)-1D-5))
+		    if (j2 < j1)
+			next
+
+		    vnew[2] = j
+		    call aclri (Memi[bufref], nc)
+		    do l = j1, j2 {
+			vold[2] = l
+			if (!pm_linenotempty (pm, vold))
+			    next
+			call pmglpi (pm, vold, Memi[bufpm], 0, nx, 0)
+			do i = 1, nc {
+			    call mw_ctrand (ctx, i-0.5D0, x1, 1)
+			    call mw_ctrand (ctx, i+0.5D0, x2, 1)
+			    i1 = max (1, nint(min(x1,x2)+1D-5))
+			    i2 = min (ncpm, nint(max(x1,x2)-1D-5))
+			    if (i2 < i1)
+				next
+			    val = Memi[bufref+i-1]
+			    do k = i1-vold[1], i2-vold[1]
+				val = max (val, Memi[bufpm+k])
+			    Memi[bufref+i-1] = val
+			}
+		    }
+		    call pmplpi (pmnew, vnew, Memi[bufref], 0, nc, PIX_SRC)
+		}
+	    }
+	    call mfree (bufref, TY_INT)
+	    call mfree (bufpm, TY_INT)
+	}
+
+	call mw_close (mw)
+	call yt_pmunmap (im)
+	im = imnew
+	call imseti (im, IM_PMDES, pmnew)
+end
+
+
+# XT_MATCH_WORLD -- Set the pixel mask to match the reference image in
+# world coordinates.  The algorithm can fail in various ways, especially
+# when higher order WCS are used.  This ideally works with images and masks
+# that are not greatly skewed in RA/DEC space.
+
+procedure yt_match_world (im, refim, maxmaskval)
+
+pointer	im			#U Pixel mask image pointer
+pointer	refim			#I Reference image pointer
+int	maxmaskval		#I Maximum mask value
+
+int	i, j, k, l, nc, nl, ncpm, nlpm, cstep, lstep, buf, nxmsi, nymsi
+int	c_im, l_im, c_ref, l_ref, c1_ref, c2_ref, l1_ref, l2_ref
+int	xmin, xmax, ymin, ymax
+double	pix_im[2], pix_ref[2], pix_tmp[2], w1[2], w2[2]
+real	x, y, icstep, ilstep, d[2], der[2,2]
+long	v[2]
+pointer	sp, bits, rl
+pointer	ba, mw_im, mw_ref, ct1, ct2, pm, xmsi, ymsi, xvec, yvec, ptr
+
+int	imstati()
+real	msieval()
+pointer	xt_baopen(), pm_open(), im_pmmapo(), imgl1i()
+pointer	mw_openim(), mw_sctran()
+bool	pm_empty()
+errchk	xt_baopen, pm_open, mw_openim, im_pmmapo, msiinit, msifit
+
+begin
+	if (im == NULL)
+	    return
+
+	# Set sizes.
+	ncpm = IM_LEN(im,1)
+	nlpm = IM_LEN(im,2)
+	nc = IM_LEN(refim,1)
+	nl = IM_LEN(refim,2)
+
+	# If the mask is empty and the sizes are the same then it does not
+	# matter if the two are actually matched in world coordinates.
+	pm = imstati (im, IM_PMDES)
+	if (pm_empty(pm) && nc == ncpm && nl == nlpm)
+	    return
+
+	# Allocate working lines.
+	call smark (sp)
+	call salloc (bits, nc, TY_INT)
+	call salloc (rl, 3+3*ncpm, TY_INT)
+
+	# Use a bit array to hold the output in memory compactly.
+	ba = xt_baopen (nc, nl, maxmaskval)
+
+	# Set logical to logical transformation through the world coordinate
+	# systems.  Use a surface fit to speed up the WCS calculations.
+
+	mw_im = mw_openim (im)
+	mw_ref = mw_openim (refim)
+
+	# First bound the reference image in world coordinates.
+	# The image is sampled and a small amount of buffer is used.
+
+	ct1 = mw_sctran (mw_ref, "logical", "world", 3)
+	cstep = 20; lstep = 20
+	icstep = (nc - 1.) / (cstep - 1.); ilstep = (nl - 1.) / (lstep - 1.)
+	w1[1] = MAX_DOUBLE; w1[2] = -MAX_DOUBLE
+	w2[1] = MAX_DOUBLE; w2[2] = -MAX_DOUBLE
+	for (pix_im[2]=1-ilstep; pix_im[2]<=nl+1+ilstep;
+	    pix_im[2]=pix_im[2]+ilstep) {
+	    for (pix_im[1]=1-icstep; pix_im[1]<=nc+1+icstep;
+		pix_im[1]=pix_im[1]+icstep) {
+		call mw_ctrand (ct1, pix_im, pix_ref, 2)
+		w1[1] = min (w1[1], pix_ref[1])
+		w1[2] = max (w1[2], pix_ref[1])
+		w2[1] = min (w2[1], pix_ref[2])
+		w2[2] = max (w2[2], pix_ref[2])
+	    }
+	}
+	call mw_ctfree (ct1)
+
+	# Fit coordinate surfaces for the mapping from the mask to the
+	# the reference image.  This is done because the WCS evaluations
+	# can be slow.  This is done on a subsample and then linear
+	# interpolation will be done.  Provide a buffer to avoid edge
+	# effects from the subsampling.  Bound the mask to what overlaps
+	# the reference image.
+
+	cstep = 10; lstep = 10; buf = 1
+
+	ct1 = mw_sctran (mw_im, "logical", "world", 3)
+	ct2 = mw_sctran (mw_ref, "world", "logical", 3)
+
+	call msiinit (xmsi, II_BILINEAR)
+	call msiinit (ymsi, II_BILINEAR)
+	nxmsi = nint ((ncpm - 1.) / cstep + 2*buf + 1)
+	nymsi = nint ((nlpm - 1.) / lstep + 2*buf + 1)
+	icstep = (nxmsi - (2.*buf + 1)) / (ncpm - 1.)
+	ilstep = (nymsi - (2.*buf + 1)) / (nlpm - 1.)
+	call malloc (xvec, nxmsi*nymsi, TY_REAL)
+	call malloc (yvec, nxmsi*nymsi, TY_REAL)
+	xmin=ncpm+1; xmax=0; ymin=nlpm+1; ymax=0
+	k = -1
+	do j = 1, nymsi {
+	    pix_im[2] = (j - (2*buf))  / ilstep + 1
+	    do i = 1, nxmsi {
+		k = k + 1
+		pix_im[1] = (i - (2*buf)) / icstep + 1
+		call mw_ctrand (ct1, pix_im, pix_tmp, 2)
+		if (pix_tmp[1] < w1[1] || pix_tmp[1] > w1[2] ||
+		    pix_tmp[2] < w2[1] || pix_tmp[2] > w2[2]) {
+		    Memr[xvec+k] = 0
+		    Memr[yvec+k] = 0
+		    next
+		}
+
+		call mw_ctrand (ct2, pix_tmp, pix_ref, 2)
+		x = pix_ref[1]
+		y = pix_ref[2]
+		if (x > 0.5 && x < nc+0.5 && y > 0.5 && y < nl+0.5) {
+		    l = max (1, min (ncpm, nint (pix_im[1])))
+		    xmin = min (xmin, l)
+		    xmax = max (xmax, l)
+		    l = max (1, min (nlpm, nint (pix_im[2])))
+		    ymin = min (ymin, l)
+		    ymax = max (ymax, l)
+		}
+
+		Memr[xvec+k] = x
+		Memr[yvec+k] = y
+	    }
+	}
+	call msifit (xmsi, Memr[xvec], nxmsi, nymsi, nxmsi) 
+	call msifit (ymsi, Memr[yvec], nxmsi, nymsi, nxmsi) 
+	call mfree (xvec, TY_REAL)
+	call mfree (yvec, TY_REAL)
+	call mw_close (mw_im)
+	call mw_close (mw_ref)
+
+	# Expand the mask bound to avoid missing the edge.
+	i = (xmin - 1) * icstep + (2*buf) - 1
+	xmin = (i - (2*buf)) / icstep + 1
+	xmin = max (1, min (ncpm, nint(xmin)))
+	i = (xmax - 1) * icstep + (2*buf) + 1.99
+	xmax = (i - (2*buf)) / icstep + 1
+	xmax = max (1, min (ncpm, nint(xmax)))
+	j = (ymin - 1) * ilstep + (2*buf) - 1
+	ymin = (j - (2*buf)) / ilstep + 1
+	ymin = max (1, min (nlpm, nint(ymin)))
+	j = (ymax - 1) * ilstep + (2*buf) + 1.99
+	ymax = (j - (2*buf)) / ilstep + 1
+	ymax = max (1, min (nlpm, nint(ymax)))
+
+	# Determine size of mask pixel in reference system.
+	# This is approximate because we don't take into account the
+	# shape of the transformed square mask pixels.
+
+	x = (xmax+xmin)/2; y = (ymax+ymin)/2
+	x = (x - 1) * icstep + (2*buf); y = (y - 1) * ilstep + (2*buf)
+	call msider (xmsi, x, y, der, 2, 2, 2)
+	d[1] = max (abs(der[2,1]), abs(der[1,2])) * icstep
+	call msider (ymsi, x, y, der, 2, 2, 2)
+	d[2] = max (abs(der[2,1]), abs(der[1,2])) * ilstep
+
+	# Go through each mask pixel and add to the new mask.
+	# This uses range lists to quickly skip good pixels.
+
+	v[1] = 1
+	do l_im = ymin, ymax {
+	    v[2] = l_im
+	    call plglri (pm, v, Memi[rl], 0, ncpm, PIX_SRC)
+	    y = (l_im - 1) * ilstep + (2*buf)
+	    ptr = rl
+	    do k = RL_FIRST, RLI_LEN(rl) {
+		ptr = ptr + RL_LENELEM
+		c_im = Memi[ptr+RL_XOFF]
+		if (c_im > xmax)
+		    next
+		if (c_im+Memi[ptr+RL_NOFF]-1 < xmin)
+		    next
+		x = (c_im - 1) * icstep + (2*buf) - icstep
+		do c_im = 1, Memi[ptr+RL_NOFF] {
+		    x = x + icstep
+		    pix_ref[1] = msieval (xmsi, x, y)
+		    pix_ref[2] = msieval (ymsi, x, y)
+		    pix_tmp[1] = max (1D0, pix_ref[1] - 0.45 * d[1])
+		    pix_tmp[2] = min (double(nc), pix_ref[1] + 0.45 * d[1])
+		    if (pix_tmp[2] < 1 || pix_tmp[1] > nc)
+		        next
+		    c1_ref = nint (pix_tmp[1])
+		    c2_ref = nint (pix_tmp[2])
+		    pix_tmp[1] = max (1D0, pix_ref[2] - 0.45 * d[2])
+		    pix_tmp[2] = min (double(nl), pix_ref[2] + 0.45 * d[2])
+		    if (pix_tmp[2] < 1 || pix_tmp[1] > nl)
+		        next
+		    l1_ref = nint (pix_tmp[1])
+		    l2_ref = nint (pix_tmp[2])
+		    do l_ref = l1_ref, l2_ref {
+			do c_ref = c1_ref, c2_ref {
+			    call xt_bapi (ba, c_ref, l_ref,
+			        Memi[ptr+RL_VOFF], 1)
+			}
+		    }
+		}
+	    }
+	}
+
+	call msifree (xmsi)
+	call msifree (ymsi)
+	call yt_pmunmap (im)
+
+	# Create a new pixel mask of the required size and populate.
+	# Do dummy image I/O to set the header.
+
+	pm = pm_open (NULL)
+	call pm_ssize (pm, 2, IM_LEN(refim,1), 27)
+	im = im_pmmapo (pm, NULL)
+	ptr = imgl1i (im)
+
+	do j = 1, nl {
+	    call xt_bagi (ba, 1, j, Memi[bits], nc)
+	    v[2] = j
+	    call pmplpi (pm, v, Memi[bits], 0, nc, PIX_SRC)
+	}
+
+	call imseti (im, IM_PMDES, pm)
+
+	call xt_baclose (ba)
+	call sfree (sp)
+end