Initial commit

20 files changed, 6853 insertions, 0 deletions

diff --git a/pkg/images/immatch/src/xregister/mkpkg b/pkg/images/immatch/src/xregister/mkpkg
new file mode 100644
index 00000000..262b721d
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/mkpkg
@@ -0,0 +1,25 @@
+# Make the XREGISTER task
+
+$checkout libpkg.a ../../../
+$update   libpkg.a
+$checkin  libpkg.a ../../../
+$exit
+
+libpkg.a:
+	rgxbckgrd.x	"xregister.h" <math/gsurfit.h>
+	rgxcolon.x	"xregister.h" <imhdr.h> <imset.h> <error.h>
+	rgxcorr.x	"xregister.h" <imhdr.h> <math/gsurfit.h> <math.h>
+	rgxdbio.x	"xregister.h"
+	rgxfft.x
+	rgxfit.x	"xregister.h" <math/iminterp.h> <mach.h> <math/nlfit.h>
+	rgxgpars.x	"xregister.h"
+	rgxicorr.x	"xregister.h" <ctype.h> <imhdr.h> <fset.h>
+	rgximshift.x	<imhdr.h> <imset.h> <math/iminterp.h>
+	rgxplot.x	<imhdr.h> <gset.h>
+	rgxppars.x	"xregister.h"
+	rgxregions.x	"xregister.h" <fset.h> <imhdr.h> <ctype.h>
+	rgxshow.x	"xregister.h"
+	rgxtools.x	"xregister.h"
+	rgxtransform.x	"xregister.h" <imhdr.h> <math.h>
+	t_xregister.x	"xregister.h" <fset.h> <gset.h> <imhdr.h> <imset.h>
+	;
diff --git a/pkg/images/immatch/src/xregister/oxregister.key b/pkg/images/immatch/src/xregister/oxregister.key
new file mode 100644
index 00000000..91064ff8
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/oxregister.key
@@ -0,0 +1,33 @@
+		Xregister Image Overlay Sub-menu
+
+
+?	Print help
+c  	Overlay the marked column of the reference image
+	with the same column of the input image
+l  	Overlay the marked line of the reference image
+	with the sname line of the input image
+x 	Overlay the marked column of the reference image
+	with the x and y lagged column of the input image
+y 	Overlay the marked line of the reference image
+	with the x and y lagged line of the input image
+v 	Overlay the marked column of the reference image
+	with the x and y shifted column of the input image
+h 	Overlay the marked line of the reference image
+	with the x and y shifted line of the input image
+q	Quit 
+
+
+		Image Overlay Sub-menu Colon Commands
+
+:c  [m] [n] 	Overlay the middle [mth] column of the reference image
+		with the mth [nth] column of the input image
+:l  [m] [n]	Overlay the middle [mth] line of the reference image
+		with the mth [nth]  line of the input image
+:x  [m] 	Overlay the middle [mth] column of the reference image
+		with the x and y lagged column of the input image
+:y  [m] 	Overlay the middle [mth] line of the reference image
+		with the x and y lagged line of the input image
+:v  [m] 	Overlay the middle [mth] column of the reference image
+		with the x and y shifted column of the input image
+:h  [m] 	Overlay the middle [mth] line of the reference image
+		with the x and y shifted line of the input image
diff --git a/pkg/images/immatch/src/xregister/rgxbckgrd.x b/pkg/images/immatch/src/xregister/rgxbckgrd.x
new file mode 100644
index 00000000..c9747ee6
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxbckgrd.x
@@ -0,0 +1,63 @@
+include <math/gsurfit.h>
+include "xregister.h"
+
+# RG_XSCALE -- Compute the background offset and x and y slope.
+
+procedure rg_xscale (xc, data, npts, nx, ny, offset, coeff)
+
+pointer	xc		#I pointer to the cross-correlation function
+real	data[ARB]	#I the input data
+int	npts		#I the number of points
+int	nx, ny		#I the dimensions of the original subraster
+real	offset		#I the input offset
+real	coeff[ARB]	#O the output coefficients
+
+int	wborder
+pointer	gs
+real	loreject, hireject, zero
+int	rg_xstati(), rg_znsum(), rg_znmedian(), rg_slope()
+real	rg_xstatr()
+
+begin
+	loreject = rg_xstatr (xc, LOREJECT)
+	hireject = rg_xstatr (xc, HIREJECT)
+	wborder = rg_xstati (xc, BORDER)
+
+	switch (rg_xstati (xc, BACKGRD)) {
+	case XC_BNONE:
+	    coeff[1] = offset
+	    coeff[2] = 0.0
+	    coeff[3] = 0.0
+	case XC_MEAN:
+	    if (rg_znsum (data, npts, zero, loreject, hireject) <= 0)
+		zero = 0.0
+	    coeff[1] = zero
+	    coeff[2] = 0.0
+	    coeff[3] = 0.0
+	case XC_MEDIAN:
+	    if (rg_znmedian (data, npts, zero, loreject, hireject) <= 0)
+		zero = 0.0
+	    coeff[1] = zero
+	    coeff[2] = 0.0
+	    coeff[3] = 0.0
+	case XC_SLOPE:
+	    call gsinit (gs, GS_POLYNOMIAL, 2, 2, GS_XNONE, 1.0, real (nx), 1.0,
+	        real (ny))
+	    if (rg_slope (gs, data, npts, nx, ny, wborder, wborder, loreject,
+	        hireject) == ERR) {
+		coeff[1] = 0.0
+		coeff[2] = 0.0
+		coeff[3] = 0.0
+	    } else {
+	        call gssave (gs, coeff)
+		coeff[1] = coeff[GS_SAVECOEFF+1]
+		coeff[2] = coeff[GS_SAVECOEFF+2]
+		coeff[3] = coeff[GS_SAVECOEFF+3]
+	    }
+	    call gsfree (gs)
+	default:
+	    coeff[1] = offset
+	    coeff[2] = 0.0
+	    coeff[3] = 0.0
+	}
+end
diff --git a/pkg/images/immatch/src/xregister/rgxcolon.x b/pkg/images/immatch/src/xregister/rgxcolon.x
new file mode 100644
index 00000000..cb007473
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxcolon.x
@@ -0,0 +1,508 @@
+include <error.h>
+include <imhdr.h>
+include <imset.h>
+include "xregister.h"
+
+# RG_XCOLON-- Procedure to process colon commands for setting the cross-
+# correlation parameters.
+
+procedure rg_xcolon (gd, xc, imr, im1, im2, db, dformat, tfd, reglist, cmdstr,
+	newdata, newcross, newcenter)
+
+pointer	gd			#I pointer to the graphics stream
+pointer	xc			#I pointer to cross-correlation structure
+pointer	imr			#I/O pointer to the reference image
+pointer	im1			#I/O pointer to the input image
+pointer	im2			#I/O pointer to the output image
+pointer	db			#I/O pointer to the shifts database file
+int	dformat			#I is the shifts file in database format 
+int	tfd			#I/O the transformations file descriptor
+pointer	reglist			#I/O pointer to the regions list
+char	cmdstr[ARB]		#I input command string
+int	newdata			#I/O new input data
+int	newcross		#I/O new cross-correlation function flag
+int	newcenter		#I/O new cross-correlation peak flag
+
+bool	streq()
+int	ncmd, creg, nreg, ival, stat
+pointer	sp, cmd, str
+real	rval
+int	strdic(), open(), nscan(), rg_xstati(), fntopnb()
+int	rg_xregions(), rg_xmkregions(), strlen()
+pointer	immap(), dtmap(), rg_xstatp()
+real	rg_xstatr()
+errchk	immap(), dtmap(), open(), fntopnb()
+
+begin
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Get the command.
+	call sscan (cmdstr)
+	call gargwrd (Memc[cmd], SZ_LINE)
+	if (Memc[cmd] == EOS) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Process the command.
+	ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, XCMDS)
+	switch (ncmd) {
+	case XCMD_REFIMAGE:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    call rg_xstats (xc, REFIMAGE, Memc[str], SZ_FNAME)
+	    if (Memc[cmd] == EOS || streq (Memc[cmd], Memc[str])) {
+		call printf ("%s: %s\n")
+		    call pargstr (KY_REFIMAGE)
+		    call pargstr (Memc[str])
+	    } else {
+		if (imr != NULL) {
+		    call imunmap (imr)
+		    imr = NULL
+		}
+		iferr {
+		    imr = immap (Memc[cmd], READ_ONLY, 0)
+		} then {
+		    call erract (EA_WARN)
+		    imr = immap (Memc[str], READ_ONLY, 0)
+		} else if (IM_NDIM(imr) > 2 || IM_NDIM(imr) != IM_NDIM(im1)) {
+		    call printf (
+		    "Image has the wrong number of dimensions\n")
+		    call imunmap (imr)
+		    imr = immap (Memc[str], READ_ONLY, 0)
+		} else {
+		    call rg_xsets (xc, REFIMAGE, Memc[cmd])
+		    newdata = YES; newcross = YES; newcenter = YES
+		}
+	    }
+
+	case XCMD_IMAGE:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    call rg_xstats (xc, IMAGE, Memc[str], SZ_FNAME)
+	    if (Memc[cmd] == EOS || streq (Memc[cmd], Memc[str])) {
+		call printf ("%s: %s\n")
+		    call pargstr (KY_IMAGE)
+		    call pargstr (Memc[str])
+	    } else {
+		if (im1 != NULL) {
+		    call imunmap (im1)
+		    im1 = NULL
+		}
+		iferr {
+		    im1 = immap (Memc[cmd], READ_ONLY, 0)
+		    call imseti (im1, IM_TYBNDRY, BT_NEAREST)
+		    if (IM_NDIM(im1) == 1)
+			call imseti (im1, IM_NBNDRYPIX, IM_LEN(im1,1))
+		    else
+		        call imseti (im1, IM_NBNDRYPIX,
+			    max (IM_LEN(im1,1), IM_LEN(im1,2)))
+		} then {
+		    call erract (EA_WARN)
+		    im1 = immap (Memc[str], READ_ONLY, 0)
+		    call imseti (im1, IM_TYBNDRY, BT_NEAREST)
+		    if (IM_NDIM(im1) == 1)
+			call imseti (im1, IM_NBNDRYPIX, IM_LEN(im1,1))
+		    else
+		        call imseti (im1, IM_NBNDRYPIX,
+			    max (IM_LEN(im1,1), IM_LEN(im1,2)))
+		} else if (IM_NDIM(im1) > 2 || IM_NDIM(im1) != IM_NDIM(imr)) {
+		    call printf (
+		        "Image has the wrong number of dimensions\n")
+		    call imunmap (im1)
+		    im1 = immap (Memc[str], READ_ONLY, 0)
+		    call imseti (im1, IM_TYBNDRY, BT_NEAREST)
+		    if (IM_NDIM(im1) == 1)
+			call imseti (im1, IM_NBNDRYPIX, IM_LEN(im1,1))
+		    else
+		        call imseti (im1, IM_NBNDRYPIX,
+			    max (IM_LEN(im1,1), IM_LEN(im1,2)))
+		} else {
+		    call rg_xsets (xc, IMAGE, Memc[cmd])
+		    newdata = YES; newcross = YES; newcenter = YES
+		}
+	    }
+
+	case XCMD_OUTIMAGE:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    call rg_xstats (xc, OUTIMAGE, Memc[str], SZ_FNAME)
+	    if (im2 == NULL || Memc[cmd] == EOS || streq (Memc[cmd],
+	        Memc[str])) {
+		call printf ("%s: %s\n")
+		    call pargstr (KY_OUTIMAGE)
+		    call pargstr (Memc[str])
+	    } else {
+		if (im2 != NULL) {
+		    call imunmap (im2)
+		    im2 = NULL
+		}
+		iferr {
+		    im2 = immap (Memc[cmd], NEW_COPY, im1)
+		} then {
+		    call erract (EA_WARN)
+		    im2 = immap (Memc[str], NEW_COPY, im1)
+		} else {
+		    call rg_xsets (xc, OUTIMAGE, Memc[cmd])
+		}
+	    }
+
+	case XCMD_DATABASE:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    call rg_xstats (xc, DATABASE, Memc[str], SZ_FNAME)
+	    if (Memc[cmd] == EOS || streq (Memc[cmd], Memc[str])) {
+		call printf ("%s:  %s\n")
+		    call pargstr (KY_DATABASE)
+		    call pargstr (Memc[str])
+	    } else {
+		if (db != NULL) {
+		    if (dformat == YES)
+		        call dtunmap (db)
+		    else
+			call close (db)
+		    db = NULL
+		}
+		iferr {
+		    if (dformat == YES)
+		        db = dtmap (Memc[cmd], APPEND)
+		    else
+			db = open (Memc[cmd], NEW_FILE, TEXT_FILE)
+		} then {
+		    call erract (EA_WARN)
+		    if (dformat == YES)
+		        db = dtmap (Memc[str], APPEND)
+		    else
+			db = open (Memc[str], APPEND, TEXT_FILE)
+		} else {
+		    call rg_xsets (xc, DATABASE, Memc[cmd])
+		}
+	    }
+
+	CASE XCMD_RECORD:
+	    call gargstr (Memc[cmd], SZ_LINE)
+	    if (Memc[cmd] == EOS) {
+		call rg_xstats (xc, RECORD, Memc[str], SZ_FNAME)
+		call printf ("%s: %s\n")
+		    call pargstr (KY_RECORD)
+		    call pargstr (Memc[str])
+	    } else 
+		call rg_xsets (xc, RECORD, Memc[cmd])
+
+	case XCMD_CREGION:
+
+	    call gargi (nreg)
+	    creg = rg_xstati (xc, CREGION)
+
+	    if (nscan() == 1 || (nreg == creg)) {
+		call printf ("%s: %d/%d")
+		    call pargstr (KY_CREGION)
+		    call pargi (creg)
+		    call pargi (rg_xstati (xc, NREGIONS))
+		call printf ("  [%d:%d,%d:%d]\n")
+		    call pargi (Memi[rg_xstatp (xc,RC1)+creg-1])
+		    call pargi (Memi[rg_xstatp (xc,RC2)+creg-1])
+		    call pargi (Memi[rg_xstatp (xc,RL1)+creg-1])
+		    call pargi (Memi[rg_xstatp (xc,RL2)+creg-1])
+
+	    } else {
+		if (nreg < 1 || nreg > rg_xstati (xc,NREGIONS)) {
+		    call printf ("Region %d is out of range\n")
+			call pargi (nreg)
+		} else {
+		    call printf (
+		        "Setting current region to %d: [%d:%d,%d:%d]\n")
+			call pargi (nreg)
+			call pargi (Memi[rg_xstatp (xc,RC1)+nreg-1])
+			call pargi (Memi[rg_xstatp (xc,RC2)+nreg-1])
+			call pargi (Memi[rg_xstatp (xc,RL1)+nreg-1])
+			call pargi (Memi[rg_xstatp (xc,RL2)+nreg-1])
+		    call rg_xseti (xc, CREGION, nreg)
+		    newdata = YES; newcross = YES; newcenter = YES
+		}
+
+	    }
+
+	case XCMD_REGIONS:
+
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    call rg_xstats (xc, REGIONS, Memc[str], SZ_FNAME)
+	    if (nscan() == 1 || streq (Memc[cmd], Memc[str]) || Memc[cmd] ==
+	        EOS) {
+		call printf ("%s [string/file]: %s\n")
+		    call pargstr (KY_REGIONS)
+		    call pargstr (Memc[str])
+	    } else {
+		if (reglist != NULL) {
+		    call fntclsb (reglist)
+		    reglist = NULL
+		}
+		iferr (reglist = fntopnb (Memc[cmd], NO))
+		    reglist = NULL
+		if (rg_xregions (reglist, imr, xc, 1) > 0) {
+		    call rg_xseti (xc, CREGION, 1)
+		    call rg_xsets (xc, REGIONS, Memc[cmd])
+		    newdata = YES; newcross = YES; newcenter = YES
+		} else {
+		    if (reglist != NULL) {
+		        call fntclsb (reglist)
+		        reglist = NULL
+		    }
+		    iferr (reglist = fntopnb (Memc[str], NO))
+		        reglist = NULL
+		    if (rg_xregions (reglist, imr, xc, 1) > 0)
+			;
+		    call rg_xsets (xc, REGIONS, Memc[str])
+		    call rg_xseti (xc, CREGION, 1)
+		}
+	    }
+
+	case XCMD_REFFILE:
+
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    call rg_xstats (xc, REFFILE, Memc[str], SZ_FNAME)
+	    if (Memc[cmd] == EOS || streq (Memc[cmd], Memc[str])) {
+		call printf ("%s:  %s\n")
+		    call pargstr (KY_REFFILE)
+		    call pargstr (Memc[str])
+	    } else {
+		if (tfd != NULL) {
+		    call close (tfd)
+		    tfd = NULL
+		}
+		iferr {
+		    tfd = open (Memc[cmd], READ_ONLY, TEXT_FILE)
+		} then {
+		    tfd = NULL
+		    call erract (EA_WARN)
+		    call rg_xsets (xc, REFFILE, "")
+		    call printf ("Coords file is undefined.\n")
+		} else
+		    call rg_xsets (xc, REFFILE, Memc[cmd])
+		newdata = YES; newcross = YES; newcenter = YES
+	    }
+
+	case XCMD_XLAG:
+	    call gargi (ival)
+	    if (nscan () ==  1) {
+		call printf ("%s = %d\n")
+		    call pargstr (KY_XLAG)
+		    call pargi (rg_xstati (xc, XLAG))
+	    } else {
+		call rg_xseti (xc, XLAG, ival)
+		newdata = YES; newcross = YES; newcenter = YES
+	    }
+
+	case XCMD_YLAG:
+	    call gargi (ival)
+	    if (nscan() ==  1) {
+		call printf ("%s = %d\n")
+		    call pargstr (KY_YLAG)
+		    call pargi (rg_xstati (xc, YLAG))
+	    } else {
+		call rg_xseti (xc, YLAG, ival)
+		newdata = YES; newcross = YES; newcenter = YES
+	    }
+
+	case XCMD_DXLAG:
+	    call gargi (ival)
+	    if (nscan() ==  1) {
+		call printf ("%s = %d\n")
+		    call pargstr (KY_DXLAG)
+		    call pargi (rg_xstati (xc, DXLAG))
+	    } else {
+		call rg_xseti (xc, DXLAG, ival)
+	    }
+
+	case XCMD_DYLAG:
+	    call gargi (ival)
+	    if (nscan() ==  1) {
+		call printf ("%s = %d\n")
+		    call pargstr (KY_DYLAG)
+		    call pargi (rg_xstati (xc, DYLAG))
+	    } else {
+		call rg_xseti (xc, DYLAG, ival)
+	    }
+
+	case XCMD_BACKGROUND:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    if (Memc[cmd] != EOS)
+	        call strcat (" ", Memc[cmd], SZ_LINE)
+	    call gargwrd (Memc[cmd+strlen(Memc[cmd])], SZ_LINE)
+	    if (Memc[cmd] == EOS) {
+		call rg_xstats (xc, BSTRING, Memc[str], SZ_FNAME)
+		call printf ("%s:  %s\n")
+		    call pargstr (KY_BACKGROUND)
+		    call pargstr (Memc[str])
+	    } else {
+		call rg_xsets (xc, BSTRING, Memc[cmd])
+		newdata = YES; newcross = YES; newcenter = YES
+	    }
+
+	case XCMD_BORDER:
+	    call gargi (ival)
+	    if (nscan() ==  1) {
+		call printf ("%s = %d\n")
+		    call pargstr (KY_BORDER)
+		    call pargi (rg_xstati (xc, BORDER))
+	    } else {
+		call rg_xseti (xc, BORDER, ival)
+		newdata = YES; newcross = YES; newcenter = YES
+	    }
+
+	case XCMD_LOREJECT:
+	    call gargr (rval)
+	    if (nscan() ==  1) {
+		call printf ("%s = %g\n")
+		    call pargstr (KY_LOREJECT)
+		    call pargr (rg_xstatr (xc, LOREJECT))
+	    } else {
+		call rg_xsetr (xc, LOREJECT, rval)
+		newdata = YES; newcross = YES; newcenter = YES
+	    }
+
+	case XCMD_HIREJECT:
+	    call gargr (rval)
+	    if (nscan() ==  1) {
+		call printf ("%s = %g\n")
+		    call pargstr (KY_HIREJECT)
+		    call pargr (rg_xstatr (xc, HIREJECT))
+	    } else {
+		call rg_xsetr (xc, HIREJECT, rval)
+		newdata = YES; newcross = YES; newcenter = YES
+	    }
+
+	case XCMD_APODIZE:
+	    call gargr (rval)
+	    if (nscan() ==  1) {
+		call printf ("%s = %g\n")
+		    call pargstr (KY_APODIZE)
+		    call pargr (rg_xstatr (xc, APODIZE))
+	    } else {
+		call rg_xsetr (xc, APODIZE, max (0.0, min (rval, 0.50)))
+		newdata = YES; newcross = YES; newcenter = YES
+	    }
+
+	case XCMD_CORRELATION:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    if (Memc[cmd] == EOS) {
+		call rg_xstats (xc, CSTRING, Memc[str], SZ_FNAME)
+		call printf ("%s = %s\n")
+		    call pargstr (KY_CORRELATION)
+		    call pargstr (Memc[str])
+	    } else {
+		stat = strdic (Memc[cmd], Memc[cmd], SZ_LINE, XC_CTYPES)
+		if (stat > 0) {
+		    call rg_xseti (xc, CFUNC, stat)
+		    call rg_xsets (xc, CSTRING, Memc[cmd])
+		    newcross = YES; newcenter = YES
+		}
+	    }
+
+	case XCMD_XWINDOW:
+	    call gargi (ival)
+	    if (nscan() ==  1) {
+		call printf ("%s = %d\n")
+		    call pargstr (KY_XWINDOW)
+		    call pargi (rg_xstati (xc, XWINDOW))
+	    } else {
+		call rg_xseti (xc, XWINDOW, ival)
+		newdata = YES; newcross = YES; newcenter = YES
+	    }
+
+	case XCMD_YWINDOW:
+	    call gargi (ival)
+	    if (nscan() ==  1) {
+		call printf ("%s = %d\n")
+		    call pargstr (KY_YWINDOW)
+		    call pargi (rg_xstati (xc, YWINDOW))
+	    } else {
+		call rg_xseti (xc, YWINDOW, ival)
+		newdata = YES; newcross = YES; newcenter = YES
+	    }
+
+	case XCMD_PEAKCENTER:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    if (Memc[cmd] == EOS) {
+		call rg_xstats (xc, PSTRING, Memc[str], SZ_FNAME)
+		call printf ("%s: %s\n")
+		    call pargstr (KY_PEAKCENTER)
+		    call pargstr (Memc[str])
+	    } else {
+		stat = strdic (Memc[cmd], Memc[cmd], SZ_LINE, XC_PTYPES)
+		if (stat > 0) {
+		    call rg_xseti (xc, PFUNC, stat)
+		    call rg_xsets (xc, PSTRING, Memc[cmd])
+		    newcenter = YES
+		}
+	    }
+
+	case XCMD_XCBOX:
+	    call gargi (ival)
+	    if (nscan() ==  1) {
+		call printf ("%s = %d\n")
+		    call pargstr (KY_XCBOX)
+		    call pargi (rg_xstati (xc, XCBOX))
+	    } else {
+		if (mod (ival, 2) == 0)
+		    ival = ival + 1
+		call rg_xseti (xc, XCBOX, ival)
+		newcenter = YES
+	    }
+
+	case XCMD_YCBOX:
+	    call gargi (ival)
+	    if (nscan() ==  1) {
+		call printf ("%s = %g\n")
+		    call pargstr (KY_YCBOX)
+		    call pargi (rg_xstati (xc, YCBOX))
+	    } else {
+		if (mod (ival, 2) == 0)
+		    ival = ival + 1
+		call rg_xseti (xc, YCBOX, ival)
+		newcenter = YES
+	    }
+
+	case XCMD_SHOW:
+	    call gdeactivate (gd, 0)
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, XSHOW)
+	    switch (ncmd) {
+	    case XSHOW_DATA:
+		call rg_xnshow (xc)
+	    case XSHOW_BACKGROUND:
+		call rg_xbshow (xc)
+	    case XSHOW_CORRELATION:
+		call rg_xxshow (xc)
+	    case XSHOW_PEAKCENTER:
+		call rg_xpshow (xc)
+	    default:
+		call rg_xshow (xc)
+	    }
+	    call greactivate (gd, 0)
+
+	case XCMD_MARK:
+	    call gdeactivate (gd, 0)
+	    if (reglist != NULL) {
+		call fntclsb (reglist)
+		reglist = NULL
+	    }
+	    if (rg_xmkregions (imr, xc, 1, MAX_NREGIONS, Memc[str],
+		SZ_LINE) <= 0) {
+		call rg_xstats (xc, REGIONS, Memc[str], SZ_LINE)
+		iferr (reglist = fntopnb (Memc[str], NO))
+		    reglist = NULL
+		if (rg_xregions (reglist, imr, xc, 1) > 0)
+		    ;
+		call rg_xsets (xc, REGIONS, Memc[str])
+		call rg_xseti (xc, CREGION, 1)
+	    } else {
+		call rg_xseti (xc, CREGION, 1)
+		call rg_xsets (xc, REGIONS, Memc[str])
+	        newdata = YES; newcross = YES; newcenter = YES
+	    }
+	    call greactivate (gd, 0)
+	default:
+	    call printf ("Unknown or ambiguous colon command\7\n")
+	}
+
+	call sfree (sp)
+end
diff --git a/pkg/images/immatch/src/xregister/rgxcorr.x b/pkg/images/immatch/src/xregister/rgxcorr.x
new file mode 100644
index 00000000..a708bf7a
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxcorr.x
@@ -0,0 +1,1034 @@
+include <imhdr.h>
+include <math.h>
+include <math/gsurfit.h>
+include "xregister.h"
+
+# RG_XCORR -- Compute the shift shift for an image relative to a reference
+# image using cross-correlation techniques.
+
+int procedure rg_xcorr (imr, im1, db, dformat, xc)
+
+pointer	imr		#I pointer to the reference image
+pointer	im1		#I pointer to the input image
+pointer	db		#I pointer to the shifts database 
+int	dformat		#I write shifts file in database format ?
+pointer	xc		#I pointer to the cross-correlation structure
+
+pointer	sp, image, imname
+real	xshift, yshift
+bool	streq()
+int	rg_xstati(), fscan(), nscan()
+errchk	rg_cross(), rg_xfile()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+	call salloc (imname, SZ_FNAME, TY_CHAR)
+	call rg_xstats (xc, IMAGE, Memc[image], SZ_FNAME)
+
+	# Initialize.
+	xshift = 0.0
+	yshift = 0.0
+
+	# Compute the average shift for the image.
+	switch (rg_xstati (xc, CFUNC)) {
+	case XC_DISCRETE, XC_DIFFERENCE, XC_FOURIER:
+
+	    # Write out the parameters.
+	    if (dformat == YES)
+	        call rg_xdbparams (db, xc)
+
+	    # Compute the cross-correlation function.
+	    call rg_cross (imr, im1, xc, NULL, xshift, yshift)
+	    call rg_xsetr (xc, TXSHIFT, xshift)
+	    call rg_xsetr (xc, TYSHIFT, yshift)
+
+	    # Write out the results for the individual regions.
+	    if (dformat == YES)
+	        call rg_xwreg (db, xc)
+
+	    # Write out the total shifts.
+	    if (dformat == YES)
+	        call rg_xdbshift (db, xc)
+	    else {
+		call fprintf (db, "%s  %g  %g\n")
+		    call pargstr (Memc[image])
+		    call pargr (xshift)
+		    call pargr (yshift)
+	    }
+
+	    # Set the x and y lags for the next picture.
+	    if (rg_xstati (xc, NREFPTS) > 0) {
+	        call rg_xseti (xc, XLAG, 0)
+	        call rg_xseti (xc, YLAG, 0)
+	    } else if (IS_INDEFI (rg_xstati (xc, DXLAG)) ||
+	        IS_INDEFI (rg_xstati (xc, DYLAG))) {
+	        call rg_xseti (xc, XLAG, nint (-xshift))
+	        call rg_xseti (xc, YLAG, nint (-yshift))
+	    } else {
+	        call rg_xseti (xc, XLAG, rg_xstati (xc, XLAG) + rg_xstati (xc,
+	            DXLAG))
+	        call rg_xseti (xc, YLAG, rg_xstati (xc, YLAG) + rg_xstati (xc,
+	            DYLAG))
+	    }
+
+	case XC_FILE:
+	    if (dformat == YES)
+	        call rg_xfile (db, xc, xshift, yshift)
+	    else {
+		if (fscan (db) != EOF) {
+		    call gargwrd (Memc[imname], SZ_FNAME)
+		    call gargr (xshift)
+		    call gargr (yshift)
+		    if (! streq (Memc[imname], Memc[image]) || nscan() != 3) {
+			xshift = 0.0
+			yshift = 0.0
+		    }
+		} else {
+		    xshift = 0.0
+		    yshift = 0.0
+		}
+	    }
+	    call rg_xsetr (xc, TXSHIFT, xshift)
+	    call rg_xsetr (xc, TYSHIFT, yshift)
+
+	default:
+	    call error (0, "The correlation function is undefined.")
+	}
+
+	call sfree (sp)
+
+	return (NO)
+end
+
+
+# RG_CROSS -- Compute the cross-correlation function for all the regions
+# using discrete, fourier, or difference techniques and compute the position
+# of its peak using one of several centering algorithms.
+
+procedure rg_cross (imr, im1, xc, gd, xavshift, yavshift)
+
+pointer	imr		#I pointer to the reference image
+pointer	im1		#I pointer to the input image
+pointer	xc		#I pointer to the cross correlation structure
+pointer	gd		#I pointer to graphics stream
+real	xavshift	#O x coord shift
+real	yavshift	#O y coord shift
+
+int	i, nregions, ngood
+pointer	pxshift, pyshift
+real	xshift, yshift
+int	rg_xstati(), rg_xcget(), rg_xfget()
+pointer	rg_xstatp()
+
+begin
+	# Get the pointers.
+	pxshift = rg_xstatp (xc, XSHIFTS)
+	pyshift = rg_xstatp (xc, YSHIFTS)
+	nregions = rg_xstati (xc, NREGIONS)
+
+	# Loop over the regions.
+	xavshift = 0.0
+	yavshift = 0.0
+	ngood = 0
+	do i = 1, nregions {
+
+	    # Compute the cross_correlation function.
+	    switch (rg_xstati (xc, CFUNC)) {
+	    case XC_DISCRETE, XC_DIFFERENCE:
+	        if (rg_xcget (xc, imr, im1, i) == ERR) {
+		    Memr[pxshift+i-1] = INDEFR
+		    Memr[pyshift+i-1] = INDEFR
+		    if (rg_xstatp (xc, XCOR) != NULL)
+		        call mfree (rg_xstatp (xc, XCOR), TY_REAL)
+		    call rg_xsetp (xc, XCOR, NULL)
+		    next
+	        }
+	    case XC_FOURIER:
+	        if (rg_xfget (xc, imr, im1, i) == ERR) {
+		    Memr[pxshift+i-1] = INDEFR
+		    Memr[pyshift+i-1] = INDEFR
+		    if (rg_xstatp (xc, XCOR) != NULL)
+		        call mfree (rg_xstatp (xc, XCOR), TY_REAL)
+		    call rg_xsetp (xc, XCOR, NULL)
+		    next
+	        }
+	    default:
+		call error (0, "The correlation function is undefined")
+	    }
+
+	    # Find the peak of the cross-correlation function.
+	    call rg_fit (xc, i, gd, xshift, yshift)
+
+	    # Accumulate the shifts.
+	    xavshift = xavshift + xshift
+	    yavshift = yavshift + yshift
+	    ngood = ngood + 1
+	}
+
+	# Compute the average shift.
+	if (ngood > 0) {
+	    xavshift = xavshift / ngood
+	    yavshift = yavshift / ngood
+	}
+end
+
+
+# RG_XFILE -- Read the average x and y shifts from the shifts database.
+
+procedure rg_xfile (db, xc, xshift, yshift)
+
+pointer	db		#I pointer to the database
+pointer	xc		#I pointer to the cross correlation structure
+real	xshift		#O shift in x
+real	yshift		#O shift in y
+
+int	rec
+pointer	sp, str
+int	dtlocate()
+real	dtgetr()
+errchk	dtlocate(), dtgetr()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	call rg_xstats (xc, RECORD, Memc[str], SZ_LINE)
+	iferr {
+	    rec = dtlocate (db, Memc[str])
+	    xshift = dtgetr (db, rec, "xshift")
+	    yshift = dtgetr (db, rec, "yshift")
+	} then {
+	    xshift = 0.0
+	    yshift = 0.0
+	}
+	
+	call sfree (sp)
+end
+
+
+# RG_ICROSS -- Compute the cross-correlation function for a given region.
+
+int procedure rg_icross (xc, imr, im1, nreg)
+
+pointer	xc		#I pointer to the cross-correlation structure
+pointer	imr		#I pointer to the reference image
+pointer	im1		#I pointer to the input image
+int	nreg		#I the index of the current region
+
+int	stat
+pointer	pxshift, pyshift
+int	rg_xstati(), rg_xcget(), rg_xfget()
+pointer	rg_xstatp()
+
+begin
+	pxshift = rg_xstatp (xc, XSHIFTS)
+	pyshift = rg_xstatp (xc, YSHIFTS)
+
+	switch (rg_xstati (xc, CFUNC)) {
+	case XC_DISCRETE, XC_DIFFERENCE:
+	    stat = rg_xcget (xc, imr, im1, nreg)
+	    if (stat == ERR) {
+		Memr[pxshift+nreg-1] = INDEFR
+		Memr[pyshift+nreg-1] = INDEFR
+		if (rg_xstatp (xc, XCOR) != NULL)
+		    call mfree (rg_xstatp (xc, XCOR), TY_REAL)
+		call rg_xsetp (xc, XCOR, NULL)
+	    }
+	case XC_FOURIER:
+	    stat = rg_xfget (xc, imr, im1, nreg)
+	    if (stat == ERR) {
+		Memr[pxshift+nreg-1] = INDEFR
+		Memr[pyshift+nreg-1] = INDEFR
+		if (rg_xstatp (xc, XCOR) != NULL)
+		    call mfree (rg_xstatp (xc, XCOR), TY_REAL)
+		call rg_xsetp (xc, XCOR, NULL)
+	    }
+	case XC_FILE:
+	    stat = OK
+	}
+
+	return (stat)
+end
+
+
+# RG_XCGET -- Compute the convolution using the discrete or difference
+# correlation functions.
+
+int procedure rg_xcget (xc, imr, im1, i)
+
+pointer	xc		#I pointer to the cross-correlation structure
+pointer	imr		#I pointer to the reference image
+pointer	im1		#I pointer to input image image
+int	i		#I index of region
+
+int	stat, xwindow, ywindow, nrimcols, nrimlines, nimcols, nimlines
+int	nrcols, nrlines, ncols, nlines
+int	xlag, ylag, nborder, rc1, rc2, rl1, rl2, c1, c2, l1, l2
+pointer	sp, str, coeff, rbuf, ibuf, xcor
+pointer	prc1, prc2, prl1, prl2, przero, prxslope, pryslope, border
+real	rxlag, rylag
+int	rg_xstati(), rg_border()
+pointer	rg_xstatp(), rg_ximget()
+real	rg_xstatr()
+
+define	nextregion_	10
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (coeff, max (GS_SAVECOEFF + 6, 9), TY_REAL)
+	rbuf = NULL
+	ibuf = NULL
+
+	# Check for regions.
+	if (i > rg_xstati (xc, NREGIONS)) {
+	    stat = ERR
+	    goto nextregion_
+	}
+
+	# Get the image sizes.
+	nrimcols = IM_LEN(imr,1)
+	if (IM_NDIM(imr) == 1)
+	    nrimlines = 1
+	else
+	    nrimlines = IM_LEN(imr,2)
+	nimcols = IM_LEN(im1,1)
+	if (IM_NDIM(im1) == 1)
+	    nimlines = 1
+	else
+	    nimlines = IM_LEN(im1,2)
+
+	# Get the reference region pointers.
+	prc1 = rg_xstatp (xc, RC1)
+	prc2 = rg_xstatp (xc, RC2)
+	prl1 = rg_xstatp (xc, RL1)
+	prl2 = rg_xstatp (xc, RL2)
+	przero = rg_xstatp (xc, RZERO)
+	prxslope = rg_xstatp (xc, RXSLOPE)
+	pryslope = rg_xstatp (xc, RYSLOPE)
+
+	# Compute the reference region limits.
+	rc1 = max (1, min (int (nrimcols), Memi[prc1+i-1]))
+	rc2 = min (int (nrimcols), max (1, Memi[prc2+i-1]))
+	rl1 = max (1, min (int (nrimlines), Memi[prl1+i-1]))
+	rl2 = min (int (nrimlines), max (1, Memi[prl2+i-1]))
+	nrcols = rc2 - rc1 + 1
+	nrlines = rl2 - rl1 + 1
+
+	# Move to the next reference region if current region is off the image.
+	if (rc1 > nrimcols || rc2 < 1 || rl1 > nrimlines || rl2 < 1) {
+	    call rg_xstats (xc, REFIMAGE, Memc[str], SZ_LINE)
+	    call eprintf (
+		"Reference section: %s[%d:%d,%d:%d] is off image.\n")
+		call pargstr (Memc[str])
+		call pargi (rc1)
+		call pargi (rc2)
+		call pargi (rl1)
+		call pargi (rl2)
+	    stat = ERR
+	    goto nextregion_
+	}
+
+	# Check the window sizes.
+	xwindow = rg_xstati (xc, XWINDOW)
+	if (nrlines == 1)
+	    ywindow = 1
+	else
+	    ywindow = rg_xstati (xc, YWINDOW)
+
+	# Move to next ref regions if current region is too small.
+	if (nrcols < xwindow || (IM_NDIM(imr) == 2 && nrlines < ywindow)) {
+	    call rg_xstats (xc, REFIMAGE, Memc[str], SZ_LINE)
+	    call eprintf (
+	    "Reference section: %s[%d:%d,%d:%d] has too few points.\n")
+		call pargstr (Memc[str])
+		call pargi (rc1)
+		call pargi (rc2)
+		call pargi (rl1)
+		call pargi (rl2)
+	    stat = ERR
+	    goto nextregion_
+	}
+
+	# Apply the transformation if defined or lag to the ref regions.
+	if (rg_xstati (xc, NREFPTS) > 0) {
+	    call rg_etransform (xc, (rc1 + rc2) / 2.0, (rl1 + rl2) / 2.0,
+		rxlag, rylag)
+	    xlag = rxlag - (rc1 + rc2) / 2.0
+	    if (ywindow == 1)
+	        ylag = 0
+	    else
+	        ylag = rylag - (rl1 + rl2) / 2.0
+	} else {
+	    xlag = rg_xstati (xc, XLAG)
+	    if (ywindow == 1)
+	        ylag = 0
+	    else
+	        ylag = rg_xstati (xc, YLAG)
+	}
+
+	# Get the input image limits.
+	c1 = rc1 + xlag - xwindow / 2
+	c2 = rc2 + xlag + xwindow / 2
+	l1 = rl1 + ylag - ywindow / 2
+	l2 = rl2 + ylag + ywindow / 2
+	ncols = c2 - c1 + 1
+	nlines = l2 - l1 + 1
+
+	# Move to the next ref region if input region is off image.
+	if (c1 > nimcols || c2 < 1 || l1 > nimlines || l2 < 1) {
+	    call rg_xstats (xc, IMAGE, Memc[str], SZ_LINE)
+	    call eprintf (
+		"Image section: %s[%d:%d,%d:%d] is off image.\n")
+		call pargstr (Memc[str])
+		call pargi (c1)
+		call pargi (c2)
+		call pargi (l1)
+		call pargi (l2)
+	    stat = ERR
+	    goto nextregion_
+	}
+
+	# Move to the next ref region if input region is less than 3 by 3.
+	if ((ncols < xwindow) || (IM_NDIM(im1) == 2 && nlines < ywindow)) {
+	    call rg_xstats (xc, IMAGE, Memc[str], SZ_LINE)
+	    call eprintf (
+		"Image section: %s[%d:%d,%d:%d] has too few points.\n")
+		call pargstr (Memc[str])
+		call pargi (c1)
+		call pargi (c2)
+		call pargi (l1)
+		call pargi (l2)
+	    stat = ERR
+	    goto nextregion_
+	}
+
+	# Get the input reference and input image data.
+	rbuf = rg_ximget (imr, rc1, rc2, rl1, rl2)
+	if (rbuf == NULL) {
+	    stat = ERR
+	    goto nextregion_
+	}
+	ibuf = rg_ximget (im1, c1, c2, l1, l2)
+	if (ibuf == NULL) {
+	    stat = ERR
+	    goto nextregion_
+	}
+
+	# Do the background subtraction.
+
+	# Compute the zero point, x slope and y slope of ref image.
+	if (IS_INDEFR(Memr[przero+i-1]) || IS_INDEFR(Memr[prxslope+i- 1]) ||
+	    IS_INDEFR(Memr[pryslope+i-1])) {
+	    if (IS_INDEFI(rg_xstati (xc, BORDER))) {
+		call rg_xscale (xc, Memr[rbuf], nrcols * nrlines, nrcols,
+		    nrlines, rg_xstatr (xc, BVALUER), Memr[coeff])
+	    } else {
+		    border = NULL
+		    nborder = rg_border (Memr[rbuf], nrcols, nrlines,
+		        max (0, nrcols - 2 * rg_xstati (xc, BORDER)),
+		        max (0, nrlines - 2 * rg_xstati (xc, BORDER)),
+			border)
+		    call rg_xscale (xc, Memr[border], nborder, nrcols,
+			nrlines, rg_xstatr (xc, BVALUER), Memr[coeff])
+		    if (border != NULL)
+		        call mfree (border, TY_REAL)
+	    }
+
+	    # Save the coefficients.
+	    Memr[przero+i-1] = Memr[coeff]
+	    Memr[prxslope+i-1] = Memr[coeff+1]
+	    Memr[pryslope+i-1] = Memr[coeff+2]
+	}
+
+	call rg_subtract (Memr[rbuf], nrcols, nrlines, Memr[przero+i-1],
+	    Memr[prxslope+i-1], Memr[pryslope+i-1])
+
+	# Compute the zero point, and the x and y slopes of input image.
+	if (IS_INDEFI(rg_xstati (xc, BORDER))) {
+	    call rg_xscale (xc, Memr[ibuf], ncols * nlines, ncols,
+		nlines, rg_xstatr (xc, BVALUE), Memr[coeff])
+	} else {
+	    border = NULL
+	    nborder = rg_border (Memr[ibuf], ncols, nlines,
+	        max (0, ncols - 2 * rg_xstati (xc, BORDER)),
+	        max (0, nlines - 2 * rg_xstati (xc, BORDER)),
+		border)
+	    call rg_xscale (xc, Memr[border], nborder, ncols, nlines,
+		rg_xstatr (xc, BVALUE), Memr[coeff])
+	    if (border != NULL)
+		call mfree (border, TY_REAL)
+	}
+
+	# Subtract the baseline.
+	call rg_subtract (Memr[ibuf], ncols, nlines, Memr[coeff],
+	    Memr[coeff+1], Memr[coeff+2])
+
+	# Apodize the data.
+	if (rg_xstatr (xc, APODIZE) > 0.0) {
+	    call rg_apodize (Memr[rbuf], nrcols, nrlines, rg_xstatr (xc,
+		APODIZE), YES)
+	    call rg_apodize (Memr[ibuf], ncols, nlines, rg_xstatr (xc,
+	        APODIZE), YES)
+	}
+
+	# Spatially filter the data with a Laplacian.
+	switch (rg_xstati (xc, FILTER)) {
+	case XC_LAPLACE:
+	    call rg_xlaplace (Memr[rbuf], nrcols, nrlines, 1.0)
+	    call rg_xlaplace (Memr[ibuf], ncols, nlines, 1.0)
+	default:
+	    ;
+	}
+
+	# Allocate space for the cross-correlation function.
+	if (rg_xstatp (xc, XCOR) == NULL) {
+	    call malloc (xcor, xwindow * ywindow, TY_REAL)
+	    call rg_xsetp (xc, XCOR, xcor)
+	} else {
+	    xcor = rg_xstatp (xc, XCOR)
+	    call realloc (xcor, xwindow * ywindow, TY_REAL)
+	    call rg_xsetp (xc, XCOR, xcor)
+	}
+
+	# Clear the correlation function.
+	call aclrr (Memr[xcor], xwindow * ywindow)
+
+	# Compute the cross-correlation function.
+	if (rg_xstati (xc, CFUNC) == XC_DISCRETE) {
+	    call rg_xconv (Memr[rbuf], nrcols, nrlines, Memr[ibuf], ncols,
+	        nlines, Memr[xcor], xwindow, ywindow)
+	} else {
+	    call rg_xdiff (Memr[rbuf], nrcols, nrlines, Memr[ibuf], ncols,
+	        nlines, Memr[xcor], xwindow, ywindow)
+	}
+
+	stat = OK
+
+nextregion_
+
+	# Free memory.
+	call sfree (sp)
+	if (rbuf != NULL)
+	    call mfree (rbuf, TY_REAL)
+	if (ibuf != NULL)
+	    call mfree (ibuf, TY_REAL)
+	if (stat == ERR)
+	    return (ERR)
+	else
+	    return (OK)
+end
+
+
+# RG_XFGET -- Compute the cross-correlation function using Fourier techniques.
+
+int procedure rg_xfget (xc, imr, im1, i)
+
+pointer	xc		#I pointer to the cross-correlation structure
+pointer	imr		#I pointer to the reference image
+pointer	im1		#I pointer to the input image
+int	i		#I index of the current region
+
+int	rc1, rc2, rl1, rl2, nrcols, nrlines, c1, c2, l1, l2, ncols, nlines
+int	nrimcols, nrimlines, nimcols, nimlines
+int	xwindow, ywindow, xlag, nxfft, nyfft, ylag, stat, nborder
+pointer	sp, str, coeff, xcor, rbuf, ibuf, fft, border
+pointer	prc1, prc2, prl1, prl2, przero, prxslope, pryslope
+real	rxlag, rylag
+int	rg_xstati(), rg_border(), rg_szfft()
+pointer	rg_xstatp(), rg_ximget()
+real	rg_xstatr()
+
+define	nextregion_	11
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (coeff, max (GS_SAVECOEFF+6, 9), TY_REAL)
+
+	# Check for number of regions.
+	if (i > rg_xstati (xc, NREGIONS)) {
+	    stat = ERR
+	    goto nextregion_
+	}
+
+	# Allocate space for the cross-correlation function.
+	nrimcols = IM_LEN(imr,1)
+	if (IM_NDIM(imr) == 1)
+	    nrimlines = 1
+	else
+	    nrimlines = IM_LEN(imr,2)
+	nimcols = IM_LEN(im1,1)
+	if (IM_NDIM(im1) == 1)
+	    nimlines = 1
+	else
+	    nimlines = IM_LEN(im1,2)
+
+	# Get the regions pointers.
+	prc1 = rg_xstatp (xc, RC1)
+	prc2 = rg_xstatp (xc, RC2)
+	prl1 = rg_xstatp (xc, RL1)
+	prl2 = rg_xstatp (xc, RL2)
+	przero = rg_xstatp (xc, RZERO)
+	prxslope = rg_xstatp (xc, RXSLOPE)
+	pryslope = rg_xstatp (xc, RYSLOPE)
+
+	# Get the reference subraster region.
+	rc1 = max (1, min (int (nrimcols), Memi[prc1+i-1]))
+	rc2 = min (int (nrimcols), max (1, Memi[prc2+i-1]))
+	rl1 = max (1, min (int (nrimlines), Memi[prl1+i-1]))
+	rl2 = min (int (nrimlines), max (1, Memi[prl2+i-1]))
+	nrcols = rc2 - rc1 + 1
+	nrlines = rl2 - rl1 + 1
+
+	# Go to next region if the reference region is off the image.
+	if (rc1 > nrimcols || rc2 < 1 || rl1 > nrimlines || rl2 < 1) {
+	    call rg_xstats (xc, REFIMAGE, Memc[str], SZ_LINE)
+	    call eprintf (
+		"Reference section: %s[%d:%d,%d:%d] is off image.\n")
+		call pargstr (Memc[str])
+		call pargi (rc1)
+		call pargi (rc2)
+		call pargi (rl1)
+		call pargi (rl2)
+	    stat = ERR
+	    goto nextregion_
+	}
+
+	# Check the window sizes.
+	xwindow = rg_xstati (xc, XWINDOW)
+	if (nrlines == 1)
+	    ywindow = 1
+	else
+	    ywindow = rg_xstati (xc, YWINDOW)
+
+	# Go to the next region if the reference region has too few points.
+	if ((nrcols < xwindow) || (IM_NDIM(im1) == 2 && nrlines < ywindow)) {
+	    call rg_xstats (xc, REFIMAGE, Memc[str], SZ_LINE)
+	    call eprintf (
+		"Reference section: %s[%d:%d,%d:%d] has too few points.\n")
+		call pargstr (Memc[str])
+		call pargi (rc1)
+		call pargi (rc2)
+		call pargi (rl1)
+		call pargi (rl2)
+	    stat = ERR
+	    goto nextregion_
+	}
+
+	# Apply the transformation if defined or the lag.
+	if (rg_xstati (xc, NREFPTS) > 0) {
+	    call rg_etransform (xc, (rc1 + rc2) / 2.0, (rl1 + rl2) / 2.0,
+		rxlag, rylag)
+	    xlag = rxlag - (rc1 + rc2) / 2.0
+	    if (ywindow == 1)
+		ylag = 0
+	    else
+	        ylag = rylag - (rl1 + rl2) / 2.0
+	} else {
+	    xlag = rg_xstati (xc, XLAG)
+	    if (ywindow == 1)
+		ylag = 0
+	    else
+	        ylag = rg_xstati (xc, YLAG)
+	}
+
+	# Get the input image subraster regions.
+	c1 = rc1 + xlag
+	c2 = rc2 + xlag
+	l1 = rl1 + ylag
+	l2 = rl2 + ylag
+	ncols = c2 - c1 + 1
+	nlines = l2 - l1 + 1
+
+	# Go to next region if region is off the image.
+	if (c1 > nimcols || c2 < 1 || l1 > nimlines || l2 < 1) {
+	    call rg_xstats (xc, IMAGE, Memc[str], SZ_LINE)
+	    call eprintf (
+		"Image section: %s[%d:%d,%d:%d] is off image.\n")
+		call pargstr (Memc[str])
+		call pargi (c1)
+		call pargi (c2)
+		call pargi (l1)
+		call pargi (l2)
+	    stat = ERR
+	    goto nextregion_
+	}
+
+	# Go to next region if region has too few points.
+	if ((ncols < xwindow) || (IM_NDIM(im1) == 2 && nlines < ywindow)) {
+	    call rg_xstats (xc, IMAGE, Memc[str], SZ_LINE)
+	    call eprintf (
+		"Image section: %s[%d:%d,%d:%d] has too few points.\n")
+		call pargstr (Memc[str])
+		call pargi (c1)
+		call pargi (c2)
+		call pargi (l1)
+		call pargi (l2)
+	    stat = ERR
+	    goto nextregion_
+	}
+
+	# Figure out how big the Fourier transform has to be, given
+	# the size of the reference subraster, the window size and
+	# the fact that the FFT must be a power of 2.
+
+	nxfft = rg_szfft (nrcols, xwindow)
+	if (ywindow == 1)
+	    nyfft = 1
+	else
+	    nyfft = rg_szfft (nrlines, ywindow)
+	call calloc (fft, 2 * nxfft * nyfft, TY_REAL)
+
+	# Get the input reference and input image data.
+	rbuf = NULL
+	rbuf = rg_ximget (imr, rc1, rc2, rl1, rl2)
+	if (rbuf == NULL) {
+	    stat = ERR
+	    goto nextregion_
+	}
+
+	# Do the background subtraction.
+
+	# Compute the zero point, x slope and y slope of ref image.
+	if (IS_INDEFR(Memr[przero+i-1]) || IS_INDEFR(Memr[prxslope+i- 1]) ||
+	    IS_INDEFR(Memr[pryslope+i-1])) {
+	    if (IS_INDEFI(rg_xstati (xc, BORDER))) {
+		call rg_xscale (xc, Memr[rbuf], nrcols * nrlines, nrcols,
+		    nrlines, rg_xstatr (xc, BVALUER), Memr[coeff])
+	    } else {
+		    border = NULL
+		    nborder = rg_border (Memr[rbuf], nrcols, nrlines,
+		        max (0, nrcols - 2 * rg_xstati (xc, BORDER)),
+		        max (0, nrlines - 2 * rg_xstati (xc, BORDER)),
+			border)
+		    call rg_xscale (xc, Memr[border], nborder, nrcols,
+			nrlines, rg_xstatr (xc, BVALUER), Memr[coeff])
+		    if (border != NULL)
+		        call mfree (border, TY_REAL)
+	    }
+
+	    # Save the coefficients.
+	    Memr[przero+i-1] = Memr[coeff]
+	    Memr[prxslope+i-1] = Memr[coeff+1]
+	    Memr[pryslope+i-1] = Memr[coeff+2]
+	}
+
+	call rg_subtract (Memr[rbuf], nrcols, nrlines, Memr[przero+i-1],
+	    Memr[prxslope+i-1], Memr[pryslope+i-1])
+
+	# Apodize the data.
+	if (rg_xstatr (xc, APODIZE) > 0.0)
+	    call rg_apodize (Memr[rbuf], nrcols, nrlines, rg_xstatr (xc,
+		APODIZE), YES)
+
+	# Spatially filter the data with a Laplacian.
+	switch (rg_xstati (xc, FILTER)) {
+	case XC_LAPLACE:
+	    call rg_xlaplace (Memr[rbuf], nrcols, nrlines, 1.0)
+	default:
+	    ;
+	}
+
+	# Load the reference data into the  FFT.
+	call rg_rload (Memr[rbuf], nrcols, nrlines, Memr[fft], nxfft, nyfft)
+	call mfree (rbuf, TY_REAL)
+
+	ibuf = NULL
+	ibuf = rg_ximget (im1, c1, c2, l1, l2)
+	if (ibuf == NULL) {
+	    stat = ERR
+	    goto nextregion_
+	}
+
+	# Compute the zero point, and the x and y slopes of input image.
+	if (IS_INDEFI(rg_xstati (xc, BORDER))) {
+	    call rg_xscale (xc, Memr[ibuf], ncols * nlines, ncols,
+		nlines, rg_xstatr (xc, BVALUE), Memr[coeff])
+	} else {
+	    border = NULL
+	    nborder = rg_border (Memr[ibuf], ncols, nlines,
+	        max (0, ncols - 2 * rg_xstati (xc, BORDER)),
+	        max (0, nlines - 2 * rg_xstati (xc, BORDER)),
+		border)
+	    call rg_xscale (xc, Memr[border], nborder, ncols, nlines,
+		rg_xstatr (xc, BVALUE), Memr[coeff])
+	    if (border != NULL)
+		call mfree (border, TY_REAL)
+	}
+
+	# Subtract the baseline.
+	call rg_subtract (Memr[ibuf], ncols, nlines, Memr[coeff],
+	    Memr[coeff+1], Memr[coeff+2])
+
+	# Apodize the data.
+	if (rg_xstatr (xc, APODIZE) > 0.0)
+	    call rg_apodize (Memr[ibuf], ncols, nlines, rg_xstatr (xc,
+	        APODIZE), YES)
+
+	# Spatially filter the data with a Laplacian.
+	switch (rg_xstati (xc, FILTER)) {
+	case XC_LAPLACE:
+	    call rg_xlaplace (Memr[ibuf], ncols, nlines, 1.0)
+	default:
+	    ;
+	}
+
+	# Load the image data into the FFT.
+	call rg_iload (Memr[ibuf], ncols, nlines, Memr[fft], nxfft, nyfft)
+	call mfree (ibuf, TY_REAL)
+
+	# Normalize the data.
+	call rg_fnorm (Memr[fft], nrcols, nrlines, nxfft, nyfft)
+
+	# Compute the cross-correlation function.
+	call rg_fftcor (Memr[fft], nxfft, nyfft)
+
+	# Allocate space for the correlation function.
+	if (rg_xstatp (xc, XCOR) == NULL) {
+	    call malloc (xcor, xwindow * ywindow, TY_REAL)
+	    call rg_xsetp (xc, XCOR, xcor)
+	} else {
+	    xcor = rg_xstatp (xc, XCOR)
+	    call realloc (xcor, xwindow * ywindow, TY_REAL)
+	    call rg_xsetp (xc, XCOR, xcor)
+	}
+
+	# Move the valid lags into the crosscorrelation array
+	call rg_movexr (Memr[fft], nxfft, nyfft, Memr[xcor], xwindow, ywindow)
+
+	# Free space.
+	call mfree (fft, TY_REAL)
+
+	stat = OK
+
+nextregion_
+
+	call sfree (sp)
+	if (stat == ERR)
+	    return (ERR)
+	else
+	    return (OK)
+end
+
+
+# RG_XIMGET -- Fill a buffer from a specified region of the image.
+
+pointer procedure rg_ximget (im, c1, c2, l1, l2)
+
+pointer	im		#I pointer to the iraf image
+int	c1, c2		#I column limits in the input image
+int	l1, l2		#I line limits in the input image
+
+int	i, ncols, nlines, npts
+pointer	ptr, index, buf
+pointer	imgs1r(), imgs2r()
+
+begin
+	ncols = c2 - c1 + 1
+	nlines = l2 - l1 + 1
+	npts = ncols * nlines
+	call malloc (ptr, npts, TY_REAL)
+
+	index = ptr
+	do i = l1, l2 {
+	    if (IM_NDIM(im) == 1)
+	        buf = imgs1r (im, c1, c2)
+	    else
+	        buf = imgs2r (im, c1, c2, i, i)
+	    call amovr (Memr[buf], Memr[index], ncols)
+	    index = index + ncols
+	}
+
+	return (ptr)
+end
+
+
+# RG_XLAPLACE -- Compute the Laplacian of an image subraster in place.
+
+procedure rg_xlaplace (data, nx, ny, rho)
+
+real	data[nx,ARB]		#I the input array
+int	nx, ny			#I the size of the input/output data array
+real	rho			#I the pixel to pixel correlation factor
+
+int	i, inline, outline, nxk, nyk, nxc
+pointer	sp, lineptrs, ptr
+real	rhosq, kernel[3,3]
+data	nxk /3/, nyk /3/
+
+begin
+	# Define the kernel.
+	rhosq = rho * rho
+	kernel[1,1] = rhosq
+	kernel[2,1] = -rho * (1.0 + rhosq)
+	kernel[3,1] = rhosq
+	kernel[1,2] = -rho * (1.0 + rhosq)
+	kernel[2,2] = (1.0 + rhosq) * (1 + rhosq)
+	kernel[3,2] = -rho * (1.0 + rhosq)
+	kernel[1,3] = rhosq
+	kernel[2,3] = -rho * (1.0 + rhosq)
+	kernel[3,3] = rhosq
+
+	# Set up an array of line pointers.
+	call smark (sp)
+	call salloc (lineptrs, nyk, TY_POINTER)
+
+	# Allocate working space.
+	nxc = nx + 2 * (nxk / 2)
+	do i = 1, nyk
+	    call salloc (Memi[lineptrs+i-1], nxc, TY_REAL) 
+
+	inline = 1 - nyk / 2
+	do i = 1, nyk - 1 {
+	    if (inline < 1) {
+	        call amovr (data[1,1], Memr[Memi[lineptrs+i]+nxk/2], nx)
+	        Memr[Memi[lineptrs+i]] = data[1,1]
+	        Memr[Memi[lineptrs+i]+nxc-1] = data[nx,1]
+	    } else {
+	        call amovr (data[1,i-1], Memr[Memi[lineptrs+i]+nxk/2], nx)
+	        Memr[Memi[lineptrs+i]] = data[1,i-1]
+	        Memr[Memi[lineptrs+i]+nxc-1] = data[nx,i-1]
+	    }
+	    inline = inline + 1
+	}
+
+	# Generate the output image line by line
+	do outline = 1, ny {
+
+	    # Scroll the input buffers
+	    ptr = Memi[lineptrs] 
+	    do i = 1, nyk - 1
+		Memi[lineptrs+i-1] = Memi[lineptrs+i]
+	    Memi[lineptrs+nyk-1] = ptr 
+
+	    # Read in new image line
+	    if (inline > ny) {
+		call amovr (data[1,ny], Memr[Memi[lineptrs+nyk-1]+nxk/2],
+		    nx)
+	        Memr[Memi[lineptrs+nyk-1]] = data[1,ny]
+	        Memr[Memi[lineptrs+nyk-1]+nxc-1] = data[nx,ny]
+	    } else {
+		call amovr (data[1,inline], Memr[Memi[lineptrs+nyk-1]+nxk/2],
+		    nx)
+	        Memr[Memi[lineptrs+nyk-1]] = data[1,inline]
+	        Memr[Memi[lineptrs+nyk-1]+nxc-1] = data[nx,inline]
+	    }
+
+	    # Generate output image line
+	    call aclrr (data[1,outline], nx)
+	    do i = 1, nyk
+		call acnvr (Memr[Memi[lineptrs+i-1]], data[1,outline], nx,
+		    kernel[1,i], nxk)
+
+	    inline = inline + 1
+	}
+
+	# Free the image buffer pointers
+	call sfree (sp)
+end
+
+
+# RG_XCONV -- Compute the cross-correlation function directly in the spatial
+# domain. 
+
+procedure rg_xconv (ref, nrcols, nrlines, image, ncols, nlines, xcor, xwindow,
+	ywindow)
+
+real	ref[nrcols,nrlines]	#I the input reference subraster
+int	nrcols, nrlines		#I size of the reference subraster
+real	image[ncols,nlines]	#I the input image subraster
+int	ncols, nlines		#I size of the image subraster
+real	xcor[xwindow,ywindow]	#O the output cross-correlation function
+int	xwindow, ywindow	#I size of the cross-correlation function
+
+int	lagx, lagy, i, j
+real	meanr, facr, meani, faci, sum
+real	asumr()
+#real	cxmin, cxmax
+
+begin
+	meanr = asumr (ref, nrcols * nrlines) / (nrcols * nrlines)
+	facr = 0.0
+	do j = 1, nrlines {
+	    do i = 1, nrcols
+		facr = facr + (ref[i,j] - meanr) ** 2
+	}
+	if (facr <= 0.0)
+	    facr = 1.0
+	else
+	    facr = sqrt (facr)
+
+	do lagy = 1, ywindow {
+	    do lagx = 1, xwindow {
+		meani = 0.0
+		do j = 1, nrlines {
+		    do i = 1, nrcols 
+			meani = meani + image[i+lagx-1,j+lagy-1]
+		}
+		meani = meani / (nrcols * nrlines)
+		faci = 0.0
+		sum = 0.0
+		do j = 1, nrlines {
+		    do i = 1, nrcols {
+			faci = faci + (image[i+lagx-1,j+lagy-1] - meani) ** 2
+			sum = sum + (ref[i,j] - meanr)  *
+			    (image[i+lagx-1,j+lagy-1] - meani)
+		    }
+		}
+		if (faci <= 0.0)
+		    faci = 1.0
+		else
+		    faci = sqrt (faci)
+		xcor[lagx,lagy] = sum / facr / faci
+	    }
+	}
+end
+
+
+# RG_XDIFF -- Compute the error function at each of several templates.
+
+procedure rg_xdiff (ref, nrcols, nrlines, image, ncols, nlines, xcor, xwindow,
+	ywindow)
+
+real	ref[nrcols,nrlines]	#I reference subraste
+int	nrcols, nrlines		#I size of the reference subraster
+real	image[ncols,nlines]	#I image subraster
+int	ncols, nlines		#I size of image subraster
+real	xcor[xwindow,ywindow]	#O crosscorrelation function
+int	xwindow, ywindow	#I size of correlation function
+
+int	lagx, lagy, i, j
+real	meanr, meani, sum, cormin, cormax
+real	asumr()
+
+
+begin
+	meanr = asumr (ref, nrcols * nrlines) / (nrcols * nrlines)
+	do lagy = 1, ywindow {
+	    do lagx = 1, xwindow {
+		meani = 0.0
+		do j = 1, nrlines {
+		    do i = 1, nrcols
+			meani = meani + image[i+lagx-1,j+lagy-1]
+		}
+		meani = meani / (nrcols * nrlines)
+		sum = 0.0
+		do j = 1, nrlines {
+		    do i = 1, nrcols {
+			sum = sum + abs ((ref[i,j] - meanr) -
+			    (image[i+lagx-1,j+lagy-1] - meani))
+		    }
+		}
+		xcor[lagx,lagy] = sum 
+	    }
+	}
+
+	call alimr (xcor, xwindow * ywindow, cormin, cormax)
+	call adivkr (xcor, cormax, xcor, xwindow * ywindow)
+	call asubkr (xcor, 1.0, xcor, xwindow * ywindow)
+	call anegr (xcor, xcor, xwindow * ywindow)
+end
+
diff --git a/pkg/images/immatch/src/xregister/rgxdbio.x b/pkg/images/immatch/src/xregister/rgxdbio.x
new file mode 100644
index 00000000..3e197636
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxdbio.x
@@ -0,0 +1,290 @@
+include "xregister.h"
+
+# RG_XWREC -- Procedure to write out the whole record.
+
+procedure rg_xwrec (db, dformat, xc)
+
+pointer	db		#I pointer to the database file
+int	dformat		#I is the shifts file in database format
+pointer	xc		#I pointer to the cross correlation structure
+
+int	i, nregions, ngood, c1, c2, l1, l2, xlag, ylag
+pointer	sp, image, prc1, prc2, prl1, prl2, pxshift, pyshift
+real	xin, yin, xout, yout, xavshift, yavshift
+int	rg_xstati()
+pointer	rg_xstatp()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+
+	# Write the header record.
+	if (dformat == YES)
+	    call rg_xdbparams (db, xc)
+
+	# Fetch the pointers to the columns.
+	prc1 = rg_xstatp (xc, RC1)
+	prc2 = rg_xstatp (xc, RC2)
+	prl1 = rg_xstatp (xc, RL1)
+	prl2 = rg_xstatp (xc, RL2)
+	pxshift = rg_xstatp (xc, XSHIFTS)
+	pyshift = rg_xstatp (xc, YSHIFTS)
+	nregions = rg_xstati (xc, NREGIONS)
+
+	xavshift = 0.0
+	yavshift = 0.0
+	ngood = 0
+	do i = 1, nregions {
+
+	    xin = (Memi[prc1+i-1] + Memi[prc2+i-1]) / 2.0
+	    yin = (Memi[prl1+i-1] + Memi[prl2+i-1]) / 2.0
+	    if (rg_xstati (xc, NREFPTS) > 0) {
+	        call rg_etransform (xc, xin, yin, xout, yout)
+		xlag = xout - xin
+		ylag = yout - yin
+	    } else {
+		xlag = rg_xstati (xc, XLAG)
+		ylag = rg_xstati (xc, YLAG)
+	    }
+	    c1 = Memi[prc1+i-1] + xlag
+	    c2 = Memi[prc2+i-1] + xlag
+	    l1 = Memi[prl1+i-1] + ylag
+	    l2 = Memi[prl2+i-1] + ylag
+
+	    if (IS_INDEFR(Memr[pxshift+i-1]) || IS_INDEFR(Memr[pyshift+i-1])) {
+		if (dformat == YES)
+	            call rg_xdbshiftr (db, Memi[prc1+i-1], Memi[prc2+i-1],
+		        Memi[prl1+i-1], Memi[prl2+i-1], c1, c2, l1, l2,
+		        INDEFR, INDEFR)
+	    } else {
+		if (dformat == YES)
+	            call rg_xdbshiftr (db, Memi[prc1+i-1], Memi[prc2+i-1],
+		        Memi[prl1+i-1], Memi[prl2+i-1], c1, c2, l1, l2,
+		        Memr[pxshift+i-1], Memr[pyshift+i-1])
+		ngood = ngood + 1
+	        xavshift = xavshift + Memr[pxshift+i-1]
+	        yavshift = yavshift + Memr[pyshift+i-1]
+	    }
+	}
+
+	# Compute the average shift.
+	if (ngood <= 0) {
+	    xavshift = 0.0
+	    yavshift = 0.0
+	} else {
+	    xavshift = xavshift / ngood
+	    yavshift = yavshift / ngood
+	}
+	call rg_xsetr (xc, TXSHIFT, xavshift)
+	call rg_xsetr (xc, TYSHIFT, yavshift)
+
+	if (dformat == YES)
+	    call rg_xdbshift (db, xc)
+	else {
+	    call rg_xstats (xc, IMAGE, Memc[image], SZ_FNAME)
+	    call fprintf (db, "%s  %g  %g\n")
+		call pargstr (Memc[image])
+		call pargr (xavshift)
+		call pargr (yavshift)
+	}
+
+	call sfree (sp)
+end
+
+
+# RG_XDBPARAMS -- Write the cross-correlation parameters to the database file.
+
+procedure rg_xdbparams (db, xc)
+
+pointer	db		#I pointer to the database file
+pointer	xc		#I pointer to the cross-correlation structure
+
+pointer	sp, str
+int	rg_xstati()
+#real	rg_xstatr()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	# Write out the time record was written.
+	call dtput (db, "\n")
+	call dtptime (db)
+
+	# Write out the record name.
+	call rg_xstats (xc, RECORD, Memc[str], SZ_FNAME)
+	call dtput (db, "begin\t%s\n")
+	    call pargstr (Memc[str])
+
+	# Write the image names.
+	call rg_xstats (xc, IMAGE, Memc[str], SZ_FNAME)
+	call dtput (db, "\t%s\t\t%s\n")
+	    call pargstr (KY_IMAGE)
+	    call pargstr (Memc[str])
+	call rg_xstats (xc, REFIMAGE, Memc[str], SZ_FNAME)
+	call dtput (db, "\t%s\t%s\n")
+	    call pargstr (KY_REFIMAGE)
+	    call pargstr (Memc[str])
+
+	call dtput (db, "\t%s\t%d\n")
+	    call pargstr (KY_NREGIONS)
+	    call pargi (rg_xstati (xc, NREGIONS))
+
+	call sfree (sp)
+end
+
+
+# RG_XWREG -- Write out the results for each region individually into
+# the shifts file.
+
+procedure rg_xwreg (db, xc)
+
+pointer	db		#I pointer to the database file
+pointer	xc		#I pointer to the cross-correlation structure
+
+int	i, nregions, c1, c2, l1, l2, xlag, ylag
+pointer	prc1, prc2, prl1, prl2, pxshift, pyshift
+real	xin, yin, xout, yout
+int	rg_xstati()
+pointer	rg_xstatp()
+
+begin
+	# Fetch the regions pointers.
+	prc1 = rg_xstatp (xc, RC1)
+	prc2 = rg_xstatp (xc, RC2)
+	prl1 = rg_xstatp (xc, RL1)
+	prl2 = rg_xstatp (xc, RL2)
+	pxshift = rg_xstatp (xc, XSHIFTS)
+	pyshift = rg_xstatp (xc, YSHIFTS)
+	nregions = rg_xstati (xc, NREGIONS)
+
+	# Write out the reference image region(s) and the equivalent
+	# input image regions.
+	do i = 1, nregions {
+
+	    xin = (Memi[prc1+i-1] + Memi[prc2+i-1]) / 2.0
+	    yin = (Memi[prl1+i-1] + Memi[prl2+i-1]) / 2.0
+	    if (rg_xstati (xc, NREFPTS) > 0) {
+	        call rg_etransform (xc, xin, yin, xout, yout)
+		xlag = xout - xin
+		ylag = yout - yin
+	    } else {
+		xlag = rg_xstati (xc, XLAG)
+		ylag = rg_xstati (xc, YLAG)
+	    }
+	    c1 = Memi[prc1+i-1] + xlag
+	    c2 = Memi[prc2+i-1] + xlag
+	    l1 = Memi[prl1+i-1] + ylag
+	    l2 = Memi[prl2+i-1] + ylag
+
+	    if (IS_INDEFR(Memr[pxshift+i-1]) || IS_INDEFR(Memr[pyshift+i-1]))
+	        call rg_xdbshiftr (db, Memi[prc1+i-1], Memi[prc2+i-1],
+		    Memi[prl1+i-1], Memi[prl2+i-1], c1, c2, l1, l2,
+		    INDEFR, INDEFR)
+	    else
+	        call rg_xdbshiftr (db, Memi[prc1+i-1], Memi[prc2+i-1],
+		    Memi[prl1+i-1], Memi[prl2+i-1], c1, c2, l1, l2,
+		    Memr[pxshift+i-1], Memr[pyshift+i-1])
+	}
+end
+
+
+# RG_XDBSHIFTR -- Write out the reference image section, input image
+# section and x and y shifts for each region.
+
+procedure rg_xdbshiftr (db, rc1, rc2, rl1, rl2, c1, c2, l1, l2, xshift, yshift)
+
+pointer	db		#I pointer to the database file
+int	rc1, rc2	#I reference region column limits
+int	rl1, rl2	#I reference region line limits
+int	c1, c2		#I image region column limits
+int	l1, l2		#I image region line limits
+real	xshift		#I x shift
+real	yshift		#I y shift
+
+begin
+	call dtput (db,"\t[%d:%d,%d:%d]\t[%d:%d,%d:%d]\t%g\t%g\n")
+	    call pargi (rc1)
+	    call pargi (rc2)
+	    call pargi (rl1)
+	    call pargi (rl2)
+	    call pargi (c1)
+	    call pargi (c2)
+	    call pargi (l1)
+	    call pargi (l2)
+	    call pargr (xshift)
+	    call pargr (yshift)
+end
+
+
+# RG_XDBSHIFT -- Write the average shifts to the shifts database.
+
+procedure rg_xdbshift (db, xc)
+
+pointer	db		#I pointer to text database file
+pointer	xc		#I pointer to the cross-correlation structure
+
+real	rg_xstatr()
+
+begin
+	call dtput (db, "\t%s\t\t%g\n")
+	    call pargstr (KY_TXSHIFT)
+	    call pargr (rg_xstatr (xc, TXSHIFT))
+	call dtput (db, "\t%s\t\t%g\n")
+	    call pargstr (KY_TYSHIFT)
+	    call pargr (rg_xstatr (xc, TYSHIFT))
+end
+
+
+# RG_XPWREC -- Print the computed shift for a region.
+
+procedure rg_xpwrec (xc, i)
+
+pointer	xc		#I pointer to the cross-correlation structure
+int	i		#I the current region
+
+int	xlag, ylag, c1, c2, l1, l2
+pointer	prc1, prc2, prl1, prl2
+real	xin, yin, rxlag, rylag
+int	rg_xstati()
+pointer	rg_xstatp()
+
+begin
+	# Fetch the pointers to the reference regions.
+	prc1 = rg_xstatp (xc, RC1)
+	prc2 = rg_xstatp (xc, RC2)
+	prl1 = rg_xstatp (xc, RL1)
+	prl2 = rg_xstatp (xc, RL2)
+
+	# Transform the reference region to the input region.
+	xin = (Memi[prc1+i-1] + Memi[prc2+i-1]) / 2.0
+	yin = (Memi[prl1+i-1] + Memi[prl2+i-1]) / 2.0
+	if (rg_xstati (xc, NREFPTS) > 0) {
+	    call rg_etransform (xc, xin, yin, rxlag, rylag)
+	    xlag = rxlag - xin
+	    ylag = rylag - yin
+	} else {
+	    xlag = rg_xstati (xc, XLAG)
+	    ylag = rg_xstati (xc, YLAG)
+	}
+
+	c1 = Memi[prc1+i-1] + xlag
+	c2 = Memi[prc2+i-1] + xlag
+	l1 = Memi[prl1+i-1] + ylag
+	l2 = Memi[prl2+i-1] + ylag
+
+	# Print the results.
+	call printf ("Region %d: [%d:%d,%d:%d]  [%d:%d,%d:%d]  %g  %g\n")
+	    call pargi (i)
+	    call pargi (Memi[prc1+i-1])
+	    call pargi (Memi[prc2+i-1])
+	    call pargi (Memi[prl1+i-1])
+	    call pargi (Memi[prl2+i-1])
+	    call pargi (c1)
+	    call pargi (c2)
+	    call pargi (l1)
+	    call pargi (l2)
+	    call pargr (Memr[rg_xstatp(xc,XSHIFTS)+i-1])
+	    call pargr (Memr[rg_xstatp(xc,YSHIFTS)+i-1])
+end
diff --git a/pkg/images/immatch/src/xregister/rgxfft.x b/pkg/images/immatch/src/xregister/rgxfft.x
new file mode 100644
index 00000000..8847cf56
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxfft.x
@@ -0,0 +1,179 @@
+# RG_FFTCOR -- Compute the FFT of the reference and image data, take their
+# product,  and compute the inverse transform to get the cross-correlation
+# function. The reference and input image are loaded into alternate memory
+# locations.
+
+procedure rg_fftcor (fft, nxfft nyfft)
+
+real	fft[ARB]	#I/O array to be  fft'd
+int	nxfft, nyfft	#I dimensions of the fft
+
+pointer	sp, dim
+
+begin
+	call smark (sp)
+	call salloc (dim, 2, TY_INT)
+
+	# Fourier transform the two arrays.
+	Memi[dim] = nxfft
+	Memi[dim+1] = nyfft
+	if (Memi[dim+1] == 1)
+	    call rg_fourn (fft, Memi[dim], 1, 1)
+	else
+	    call rg_fourn (fft, Memi[dim], 2, 1)
+
+	# Compute the product of the two transforms.
+	call rg_mulfft (fft, fft, 2 * nxfft, nyfft)
+
+	# Shift the array to center the transform.
+	call rg_fshift (fft, fft, 2 * nxfft, nyfft)
+
+	# Normalize the transform.
+	call adivkr (fft, real (nxfft * nyfft), fft, 2 * nxfft * nyfft)
+
+	# Compute the inverse transform.
+	if (Memi[dim+1] == 1)
+	    call rg_fourn (fft, Memi[dim], 1, -1)
+	else
+	    call rg_fourn (fft, Memi[dim], 2, -1)
+
+	call sfree (sp)
+end
+
+
+# RG_MULFFT -- Unpack the two individual ffts and compute their product.
+
+procedure rg_mulfft (fft1, fft2, nxfft, nyfft)
+
+real	fft1[nxfft,nyfft]	#I array containing 2 ffts of 2 real functions
+real	fft2[nxfft,nyfft]	#O fft of correlation function
+int	nxfft, nyfft		#I dimensions of fft
+
+int	i,j, nxd2p2, nxp2, nxp3, nyd2p1, nyp2
+real	c1, c2, h1r, h1i, h2r, h2i
+
+begin
+	c1 = 0.5
+	c2 = -0.5
+
+	nxd2p2 = nxfft / 2 + 2
+	nxp2 = nxfft + 2
+	nxp3 = nxfft + 3
+	nyd2p1 = nyfft / 2 + 1
+	nyp2 = nyfft + 2
+
+	# Compute the 0 frequency point.
+	h1r = fft1[1,1]
+	h1i = 0.0
+	h2r = fft1[2,1]
+	h2i = 0.0
+	fft2[1,1] = h1r * h2r
+	fft2[2,1] = 0.0
+
+	# Compute the x axis points.
+	do i = 3, nxd2p2, 2 {
+	    h2r = c1 * (fft1[i,1] + fft1[nxp2-i,1])
+	    h2i = c1 * (fft1[i+1,1] - fft1[nxp3-i,1])
+	    h1r = -c2 * (fft1[i+1,1] + fft1[nxp3-i,1])
+	    h1i = c2 * (fft1[i,1] - fft1[nxp2-i,1])
+	    fft2[i,1] = (h1r * h2r + h1i * h2i)
+	    fft2[i+1,1] = (h1i * h2r - h2i * h1r)
+	    fft2[nxp2-i,1] = fft2[i,1]
+	    fft2[nxp3-i,1] = - fft2[i+1,1]
+	}
+
+	# Quit if the transform is 1D.
+	if (nyfft < 2)
+	    return
+
+	# Compute the y axis points.
+	do i = 2, nyd2p1 {
+	    h2r = c1 * (fft1[1,i] + fft1[1, nyp2-i])
+	    h2i = c1 * (fft1[2,i] - fft1[2,nyp2-i])
+	    h1r = -c2 * (fft1[2,i] + fft1[2,nyp2-i])
+	    h1i = c2 * (fft1[1,i] - fft1[1,nyp2-i])
+	    fft2[1,i] = (h1r * h2r + h1i * h2i)
+	    fft2[2,i] = (h1i * h2r - h2i * h1r)
+	    fft2[1,nyp2-i] = fft2[1,i]
+	    fft2[2,nyp2-i] = - fft2[2,i]
+	}
+
+	# Compute along the axis of symmetry.
+	do i = 3, nxd2p2, 2 {
+	    h2r = c1 * (fft1[i,nyd2p1] + fft1[nxp2-i, nyd2p1])
+	    h2i = c1 * (fft1[i+1,nyd2p1] - fft1[nxp3-i,nyd2p1])
+	    h1r = -c2 * (fft1[i+1,nyd2p1] + fft1[nxp3-i,nyd2p1])
+	    h1i = c2 * (fft1[i,nyd2p1] - fft1[nxp2-i,nyd2p1])
+	    fft2[i,nyd2p1] = (h1r * h2r + h1i * h2i)
+	    fft2[i+1,nyd2p1] = (h1i * h2r - h2i * h1r)
+	    fft2[nxp2-i,nyd2p1] = fft2[i,nyd2p1]
+	    fft2[nxp3-i,nyd2p1] = - fft2[i+1,nyd2p1]
+	}
+
+	# Compute the remainder of the transform.
+	do j = 2, nyd2p1 - 1 {
+	    do i = 3, nxfft, 2 {
+	        h2r = c1 * (fft1[i,j] + fft1[nxp2-i, nyp2-j])
+	        h2i = c1 * (fft1[i+1,j] - fft1[nxp3-i,nyp2-j])
+	        h1r = -c2 * (fft1[i+1,j] + fft1[nxp3-i,nyp2-j])
+	        h1i = c2 * (fft1[i,j] - fft1[nxp2-i,nyp2-j])
+	        fft2[i,j] = (h1r * h2r + h1i * h2i)
+	        fft2[i+1,j] = (h1i * h2r - h2i * h1r)
+	        fft2[nxp2-i,nyp2-j] = fft2[i,j]
+	        fft2[nxp3-i,nyp2-j] = - fft2[i+1,j]
+	    }
+	}
+end
+
+
+# RG_FNORM -- Normalize the reference and image data before computing
+# the fft's.
+
+procedure rg_fnorm (array, ncols, nlines, nxfft, nyfft)
+
+real	array[ARB]			#I/O the input/output data array
+int	ncols, nlines			#I dimensions of the input data array
+int	nxfft, nyfft			#I dimensions of the fft
+
+int	i, j, index
+real	sumr, sumi, meanr, meani
+
+begin
+	# Compute the mean.
+	sumr = 0.0
+	sumi = 0.0
+	index = 0
+	do j = 1, nlines {
+	    do i = 1, ncols {
+		sumr = sumr + array[index+2*i-1]
+		sumi = sumi + array[index+2*i]
+	    }
+	    index = index + 2 * nxfft
+	}
+	meanr = sumr / (ncols * nlines)
+	meani = sumi / (ncols * nlines)
+
+	# Compute the sigma.
+	sumr = 0.0
+	sumi = 0.0
+	index = 0
+	do j = 1, nlines {
+	    do i = 1, ncols {
+		sumr = sumr + (array[index+2*i-1] - meanr) ** 2
+		sumi = sumi + (array[index+2*i] - meani) ** 2
+	    }
+	    index = index + 2 * nxfft
+	}
+	sumr = sqrt (sumr)
+	sumi = sqrt (sumi)
+
+	# Normalize the data.
+	index = 0
+	do j = 1, nlines {
+	    do i = 1, ncols {
+		array[index+2*i-1] = (array[index+2*i-1] - meanr) / sumr
+		array[index+2*i] = (array[index+2*i] - meani) / sumi
+	    }
+	    index = index + 2 * nxfft
+	}
+end
diff --git a/pkg/images/immatch/src/xregister/rgxfit.x b/pkg/images/immatch/src/xregister/rgxfit.x
new file mode 100644
index 00000000..34e6398c
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxfit.x
@@ -0,0 +1,814 @@
+include <mach.h>
+include <math/iminterp.h>
+include <math/nlfit.h>
+include "xregister.h"
+
+define	NL_MAXITER		10
+define	NL_TOL			0.001
+
+# RG_FIT -- Fit the peak of the cross-correlation function using one of the
+# fitting functions.
+
+procedure rg_fit (xc, nreg, gd, xshift, yshift)
+
+pointer	xc		    #I the pointer to the cross-corrrelation structure
+int	nreg		    #I the current region
+pointer	gd		    #I the pointer to the graphics stream
+real	xshift, yshift	    #O the computed shifts
+
+int	nrlines, xwindow, ywindow, xcbox, ycbox, xlag, ylag
+real	xin, yin, xout, yout
+int	rg_xstati()
+pointer	rg_xstatp()
+
+begin
+	# Check the window and centering box sizes.
+	nrlines = Memi[rg_xstatp(xc,RL2)+nreg-1] -
+	    Memi[rg_xstatp(xc,RL1)+nreg-1] + 1
+	xwindow = rg_xstati (xc, XWINDOW)
+	if (nrlines == 1)
+	    ywindow = 1
+	else
+	    ywindow = rg_xstati (xc, YWINDOW)
+	xcbox = rg_xstati (xc, XCBOX)
+	if (nrlines == 1)
+	    ycbox = 1
+	else
+	    ycbox = rg_xstati (xc, YCBOX)
+
+	# Do the centering.
+	switch (rg_xstati (xc, PFUNC)) {
+	case XC_PNONE:
+	    call rg_maxmin (Memr[rg_xstatp(xc,XCOR)], xwindow, ywindow,
+	        xshift, yshift)
+	case XC_CENTROID:
+	    call rg_imean (Memr[rg_xstatp(xc,XCOR)], xwindow,
+	        ywindow, xcbox, ycbox, xshift, yshift)
+	case XC_SAWTOOTH:
+	    call rg_sawtooth (Memr[rg_xstatp(xc,XCOR)], xwindow,
+	        ywindow, xcbox, ycbox, xshift, yshift)
+	case XC_PARABOLA:
+	    call rg_iparabolic (Memr[rg_xstatp(xc,XCOR)], xwindow, ywindow, 
+		xcbox, ycbox, xshift, yshift)
+	case XC_MARK:
+	    if (gd == NULL)
+	        call rg_imean (Memr[rg_xstatp(xc,XCOR)], xwindow,
+	            ywindow, xcbox, ycbox, xshift, yshift)
+	    else
+	        call rg_xmkpeak (gd, xwindow, ywindow, xshift, yshift)
+	default:
+	    call rg_imean (Memr[rg_xstatp(xc,XCOR)], xwindow, ywindow,
+	        xcbox, ycbox, xshift, yshift)
+	}
+
+	# Store the shifts.
+	if (rg_xstati (xc, NREFPTS) > 0) {
+	    xin  = (Memi[rg_xstatp(xc,RC1)+nreg-1] +
+	        Memi[rg_xstatp(xc,RC2)+nreg-1]) / 2.0
+	    yin  = (Memi[rg_xstatp(xc,RL1)+nreg-1] +
+	        Memi[rg_xstatp(xc,RL2)+nreg-1]) / 2.0
+	    call rg_etransform (xc, xin, yin, xout, yout)
+	    xlag = xout - xin
+	    ylag = yout - yin
+	} else {
+	    xlag = rg_xstati (xc, XLAG)
+	    ylag = rg_xstati (xc, YLAG)
+	}
+	xshift = - (xshift + xlag)
+	yshift = - (yshift + ylag)
+	Memr[rg_xstatp(xc,XSHIFTS)+nreg-1] = xshift
+	Memr[rg_xstatp(xc,YSHIFTS)+nreg-1] = yshift
+end
+
+
+# RG_MAXMIN -- Procedure to compute the peak of the cross-correlation function
+# by determining the maximum point.
+
+procedure rg_maxmin (xcor, xwindow, ywindow, xshift, yshift)
+
+real	xcor[xwindow,ywindow]	#I the cross-correlation function
+int	xwindow, ywindow	#I dimensions of cross-correlation function
+real	xshift, yshift		#O x and shift of the peak
+
+int	xindex, yindex
+
+begin
+	# Locate the maximum point.
+	call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex)
+	xshift = xindex - (1.0 + xwindow) / 2.0
+	yshift = yindex - (1.0 + ywindow) / 2.0
+end
+
+
+# RG_IMEAN -- Compute the peak of the cross-correlation function using the
+# intensity weighted mean of the marginal distributions in x and y.
+
+procedure rg_imean (xcor, xwindow, ywindow, xcbox, ycbox, xshift, yshift)
+
+real	xcor[xwindow,ARB]	#I the cross-correlation function
+int	xwindow, ywindow	#I dimensions of the cross-correlation function
+int	xcbox, ycbox		#I dimensions of the centering box
+real	xshift, yshift		#O x and y shift of cross-correlation function
+
+int	xindex, yindex, xlo, xhi, ylo, yhi, nx, ny
+pointer	sp, xmarg, ymarg
+
+begin
+	call smark (sp)
+	call salloc (xmarg, xcbox, TY_REAL)
+	call salloc (ymarg, ycbox, TY_REAL)
+
+	# Locate the maximum point and normalize.
+	call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex)
+
+	# Compute the limits of the centering box.
+	xlo = max (1, xindex - xcbox / 2)
+	xhi = min (xwindow, xindex + xcbox / 2)
+	nx = xhi - xlo + 1
+	ylo = max (1, yindex - ycbox / 2)
+	yhi = min (ywindow, yindex + ycbox / 2)
+	ny = yhi - ylo + 1
+
+	# Accumulate the marginals.
+	call rg_xmkmarg (xcor, xwindow, ywindow, xlo, xhi, ylo, yhi,
+	    Memr[xmarg], Memr[ymarg])
+
+	# Compute the shifts.
+	call rg_centroid (Memr[xmarg], nx, xshift)
+	xshift = xshift + xlo - 1 - (1.0 + xwindow) / 2.0
+	call rg_centroid (Memr[ymarg], ny, yshift)
+	yshift = yshift + ylo - 1 - (1.0 + ywindow) / 2.0
+
+	call sfree (sp)
+end
+
+
+# RG_IPARABOLIC -- Computer the peak of the cross-correlation function by
+# doing parabolic interpolation around the peak.
+
+procedure rg_iparabolic (xcor, xwindow, ywindow, xcbox, ycbox, xshift, yshift)
+
+real	xcor[xwindow,ARB]	#I the cross-correlation function
+int	xwindow, ywindow	#I dimensions of the cross-correlation fucntion
+int	xcbox, ycbox		#I the dimensions of the centering box
+real	xshift, yshift		#O the x and y shift of the peak
+
+int	i, j, xindex, yindex, xlo, xhi, nx, ylo, yhi, ny
+pointer	sp, x, y, c, xfit, yfit
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (x, 3, TY_REAL)
+	call salloc (y, 3, TY_REAL)
+	call salloc (c, 3, TY_REAL)
+	call salloc (xfit, 3, TY_REAL)
+	call salloc (yfit, 3, TY_REAL)
+
+	# Locate the maximum point.
+	call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex)
+
+	xlo = max (1, xindex - 1)
+	xhi = min (xwindow, xindex + 1)
+	nx = xhi - xlo + 1
+	ylo = max (1, yindex - 1)
+	yhi = min (ywindow, yindex + 1)
+	ny = yhi - ylo + 1
+
+	 # Initialize.
+	 do i = 1, 3
+	    Memr[x+i-1] = i
+
+	# Fit the x shift.
+	if (nx >= 3) {
+	    do j = ylo, yhi {
+		 do i = xlo, xhi
+		    Memr[y+i-xlo] = xcor[i,j]
+		 call rg_iparab (Memr[x], Memr[y], Memr[c])
+		 Memr[xfit+j-ylo] = - Memr[c+1] / (2.0 * Memr[c+2])
+		 Memr[yfit+j-ylo] = Memr[c] + Memr[c+1] * Memr[xfit+j-ylo] +
+		     Memr[c+2] * Memr[xfit+j-ylo] ** 2 
+	    }
+	    if (ny >= 3)
+	        call rg_iparab (Memr[xfit], Memr[yfit], Memr[c])
+	    xshift = - Memr[c+1] / (2.0 * Memr[c+2])
+	} else 
+	    xshift = xindex - xlo + 1 
+
+	# Fit the y shift.
+	if (ny >= 3) {
+	    do i = xlo, xhi {
+		do j = ylo, yhi
+		    Memr[y+j-ylo] = xcor[i,j]
+		call rg_iparab (Memr[x], Memr[y], Memr[c])
+		Memr[xfit+i-xlo] = - Memr[c+1] / (2.0 * Memr[c+2])
+		Memr[yfit+i-xlo] = Memr[c] + Memr[c+1] * Memr[xfit+i-xlo] +
+		    Memr[c+2] * Memr[xfit+i-xlo] ** 2 
+	    }
+	    call rg_iparab (Memr[xfit], Memr[yfit], Memr[c])
+	    yshift = - Memr[c+1] / (2.0 * Memr[c+2])
+	} else 
+	    yshift = yindex - ylo + 1 
+
+	xshift = xshift + xlo - 1 - (1.0 + xwindow) / 2.0
+	yshift = yshift + ylo - 1 - (1.0 + ywindow) / 2.0
+
+	call sfree (sp)
+end
+
+
+define	NPARS_PARABOLA	3
+
+# RG_PARABOLIC -- Compute the peak of the cross-correlation function by fitting
+# a parabola to the peak.
+
+procedure rg_parabolic (xcor, xwindow, ywindow, xcbox, ycbox, xshift, yshift)
+
+real	xcor[xwindow,ARB]	#I the cross-correlation function
+int	xwindow, ywindow	#I dimensions of the cross-correlation fucntion
+int	xcbox, ycbox		#I the dimensions of the centering box
+real	xshift, yshift		#O the x and y shift of the peak
+
+extern	rg_polyfit, rg_dpolyfit
+int	i,  xindex, yindex, xlo, xhi, ylo, yhi, nx, ny, npar, ier
+pointer	sp, x, w, xmarg, ymarg, params, eparams, list, nl
+int	locpr()
+
+begin
+	call smark (sp)
+	call salloc (x, max (xwindow, ywindow), TY_REAL)
+	call salloc (w, max (xwindow, ywindow), TY_REAL)
+	call salloc (xmarg, max (xwindow, ywindow), TY_REAL)
+	call salloc (ymarg, max (xwindow, ywindow), TY_REAL)
+	call salloc (params, NPARS_PARABOLA, TY_REAL)
+	call salloc (eparams, NPARS_PARABOLA, TY_REAL)
+	call salloc (list, NPARS_PARABOLA, TY_INT)
+
+	# Locate the maximum point.
+	call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex)
+
+	xlo = max (1, xindex - xcbox / 2)
+	xhi = min (xwindow, xindex + xcbox / 2)
+	nx = xhi - xlo + 1
+	ylo = max (1, yindex - ycbox / 2)
+	yhi = min (ywindow, yindex + ycbox / 2)
+	ny = yhi - ylo + 1
+
+	# Accumulate the marginals.
+	call rg_xmkmarg (xcor, xwindow, ywindow, xlo, xhi, ylo, yhi,
+	    Memr[xmarg], Memr[ymarg])
+
+	# Compute the x shift.
+	if (nx >= 3) {
+	    do i = 1, nx
+		Memr[x+i-1] = i
+	    do i = 1, nx
+		Memr[w+i-1] = Memr[xmarg+i-1]
+	    call rg_iparab (Memr[x+xindex-xlo-1], Memr[xmarg+xindex-xlo-1],
+		Memr[params])
+		    xshift = - Memr[params+1] / (2.0 * Memr[params+2]) 
+		    call eprintf ("\txshift=%g\n")
+			call pargr (xshift)
+	    call aclrr (Memr[eparams], NPARS_PARABOLA)
+	    do i = 1, NPARS_PARABOLA
+		Memi[list+i-1] = i
+	    call nlinitr (nl, locpr (rg_polyfit), locpr (rg_dpolyfit),
+	        Memr[params], Memr[eparams], NPARS_PARABOLA, Memi[list],
+		    NPARS_PARABOLA, .0001, NL_MAXITER)
+	    call nlfitr (nl, Memr[x], Memr[xmarg], Memr[w], nx, 1, WTS_USER,
+		ier)
+	    call nlvectorr (nl, Memr[x], Memr[w], nx, 1)
+	    do i = 1, nx {
+		call eprintf ("x=%g y=%g yfit=%g\n")
+		    call pargr (Memr[x+i-1])
+		    call pargr (Memr[xmarg+i-1])
+		    call pargr (Memr[w+i-1])
+	    }
+	    if (ier != NO_DEG_FREEDOM) {
+		call nlpgetr (nl, Memr[params], npar)
+		if (Memr[params+2] != 0)
+		    xshift = - Memr[params+1] / (2.0 * Memr[params+2])
+		else
+	            xshift = xindex - xlo + 1 
+	    } else
+	        xshift = xindex - xlo + 1 
+	    call nlfreer (nl)
+	} else
+	    xshift = xindex - xlo + 1 
+
+	# Compute the y shift.
+	if (ny >= 3) {
+	    do i = 1, ny
+		Memr[x+i-1] = i
+	    do i = 1, ny
+		Memr[w+i-1] = Memr[ymarg+i-1]
+	    call rg_iparab (Memr[x+yindex-ylo-1], Memr[ymarg+yindex-ylo-1],
+		Memr[params])
+		    yshift = - Memr[params+1] / (2.0 * Memr[params+2]) 
+		    call eprintf ("\tyshift=%g\n")
+			call pargr (yshift)
+	    call aclrr (Memr[eparams], NPARS_PARABOLA)
+	    do i = 1, NPARS_PARABOLA
+		Memi[list+i-1] = i
+	    call nlinitr (nl, locpr (rg_polyfit), locpr (rg_dpolyfit),
+	        Memr[params], Memr[eparams], NPARS_PARABOLA, Memi[list],
+		    NPARS_PARABOLA, 0.0001, NL_MAXITER)
+	    call nlfitr (nl, Memr[x], Memr[ymarg], Memr[w], ny, 1, WTS_USER,
+		ier)
+	    call nlvectorr (nl, Memr[x], Memr[w], ny, 1)
+	    do i = 1, ny {
+		call eprintf ("x=%g y=%g yfit=%g\n")
+		    call pargr (Memr[x+i-1])
+		    call pargr (Memr[ymarg+i-1])
+		    call pargr (Memr[w+i-1])
+	    }
+	    if (ier != NO_DEG_FREEDOM) {
+		call nlpgetr (nl, Memr[params], npar)
+		if (Memr[params+2] != 0)
+		    yshift = -Memr[params+1] / (2.0 * Memr[params+2])
+		else
+	            yshift = yindex - ylo + 1 
+	    } else
+	        yshift = yindex - ylo + 1 
+	    call nlfreer (nl)
+	} else
+	    yshift = yindex - ylo + 1 
+
+	xshift = xshift + xlo - 1 - (1.0 + xwindow) / 2.0
+	yshift = yshift + ylo - 1 - (1.0 + ywindow) / 2.0
+
+	call sfree (sp)
+end
+
+define	EMISSION	1		# emission features
+define	ABSORPTION	2		# emission features
+
+# RG_SAWTOOTH -- Compute the the x and y centers using a sawtooth
+# convolution function.
+
+procedure rg_sawtooth (xcor, xwindow, ywindow, xcbox, ycbox, xshift, yshift)
+
+real	xcor[xwindow,ARB]	#I the cross-correlation function
+int	xwindow, ywindow	#I the dimensions of the cross-correlation
+int	xcbox, ycbox		#I the dimensions of the centering box
+real	xshift, yshift		#O the x and y shifts
+
+int	i, j, xindex, yindex, xlo, xhi, ylo, yhi, nx, ny
+pointer	sp, data, xfit, yfit, yclean
+real	ic
+
+begin
+	call smark (sp)
+	call salloc (data, max (xwindow, ywindow), TY_REAL)
+	call salloc (xfit, max (xwindow, ywindow), TY_REAL)
+	call salloc (yfit, max (xwindow, ywindow), TY_REAL)
+	call salloc (yclean, max (xwindow, ywindow), TY_REAL)
+
+	# Locate the maximum point and normalize.
+	call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex)
+
+	xlo = max (1, xindex - xcbox)
+	xhi = min (xwindow, xindex + xcbox)
+	nx = xhi - xlo + 1
+	ylo = max (1, yindex - ycbox)
+	yhi = min (ywindow, yindex + ycbox)
+	ny = yhi - ylo + 1
+
+	# Compute the y shift.
+	if (ny >= 3) {
+	    do j = ylo, yhi {
+	        do i = xlo, xhi
+		    Memr[data+i-xlo] = xcor[i,j]
+	        call rg_x1dcenter (real (xindex - xlo + 1), Memr[data], nx,
+	            Memr[xfit+j-ylo], Memr[yfit+j-ylo], real (nx / 2.0),
+		    EMISSION, real (nx / 2.0), 0.0)
+	    }
+	    call arbpix (Memr[yfit], Memr[yclean], ny, II_SPLINE3,
+	        II_BOUNDARYEXT) 
+	    call rg_x1dcenter (real (yindex - ylo + 1), Memr[yclean], ny,
+	        yshift, ic, real (ny / 2.0), EMISSION, real (ny / 2.0), 0.0)
+	    if (IS_INDEFR(yshift))
+	        yshift = yindex - ylo + 1
+	} else 
+	    yshift = yindex - ylo + 1
+	yshift = yshift + ylo - 1 - (1.0 + ywindow) / 2.0
+
+	# Compute the x shift.
+	if (nx >= 3) {
+	    if (ny >= 3) {
+	        do i = xlo, xhi {
+	            do j = ylo, yhi
+		        Memr[data+j-ylo] = xcor[i,j]
+	            call rg_x1dcenter (real (yindex - ylo + 1), Memr[data], ny,
+	                Memr[xfit+i-xlo], Memr[yfit+i-xlo], real (ny / 2.0),
+		        EMISSION, real (ny / 2.0), 0.0)
+	        }
+	        call arbpix (Memr[yfit], Memr[yclean], nx, II_SPLINE3,
+	            II_BOUNDARYEXT) 
+	        call rg_x1dcenter (real (xindex - xlo + 1), Memr[yclean], nx,
+	            xshift, ic, real (nx / 2.0), EMISSION, real (nx / 2.0), 0.0)
+	    } else {
+		call rg_x1dcenter (real (xindex - xlo + 1), xcor[xlo,1], nx,
+		    xshift, ic, real (nx / 2.0), EMISSION, real (nx / 2.0), 0.0)
+	    }
+	    if (IS_INDEFR(xshift))
+	        xshift = xindex - xlo + 1
+	} else
+	    xshift = xindex - xlo + 1
+	xshift = xshift + xlo - 1 - (1.0 + xwindow) / 2.0
+
+	call sfree (sp)
+end
+
+
+# RG_ALIM2R -- Determine the pixel position of the data maximum.
+
+procedure rg_alim2r (data, nx, ny, i, j)
+
+real	data[nx,ARB]		#I the input data
+int	nx, ny			#I the dimensions of the input array
+int	i, j			#O the indices of the maximum pixel
+
+int	ii, jj
+real	datamax
+
+begin
+	datamax = -MAX_REAL
+	do jj = 1, ny {
+	    do ii = 1, nx {
+		if (data[ii,jj] > datamax) {
+		    datamax = data[ii,jj]
+		    i = ii
+		    j = jj
+		}
+	    }
+	}
+end
+
+
+# RG_XMKMARG -- Acumulate the marginal arrays in x and y.
+
+procedure rg_xmkmarg (xcor, xwindow, ywindow, xlo, xhi, ylo, yhi, xmarg,
+	ymarg)
+
+real	xcor[xwindow,ARB]	#I the cross-correlation function
+int	xwindow, ywindow	#I dimensions of cross-correlation function
+int	xlo, xhi		#I the x limits for centering
+int	ylo, yhi		#I the y limits for centering
+real	xmarg[ARB]		#O the output x marginal array
+real	ymarg[ARB]		#O the output y marginal array
+
+int	i, j, index, nx, ny
+
+begin
+	nx = xhi - xlo + 1
+	ny = yhi - ylo + 1
+
+	# Compute the x marginal.
+	index = 1 - xlo
+	do i = xlo, xhi {
+	    xmarg[index+i] = 0.0
+	    do j = ylo, yhi
+		xmarg[index+i] = xmarg[index+i] + xcor[i,j]
+	}
+
+	# Normalize the x marginal.
+	call adivkr (xmarg, real (ny), xmarg, nx)
+
+	# Compute the y marginal.
+	index = 1 - ylo
+	do j = ylo, yhi {
+	    ymarg[index+j] = 0.0
+	    do i = xlo, xhi
+		ymarg[index+j] = ymarg[index+j] + xcor[i,j]
+	}
+
+	# Normalize the ymarginal.
+	call adivkr (ymarg, real (nx), ymarg, ny)
+end
+
+
+# RG_CENTROID -- Compute the intensity weighted maximum of an array.
+
+procedure rg_centroid (a, npts, shift)
+
+real	a[ARB]		#I the input array
+int	npts		#I the number of points
+real	shift		#O the position of the maximum
+
+int	i
+real	mean, dif, sumi, sumix
+bool	fp_equalr()
+real	asumr()
+
+begin
+	sumi = 0.0
+	sumix = 0.0
+	mean = asumr (a, npts) / npts
+
+	do i = 1, npts {
+	    dif = a[i]
+	    dif = a[i] - mean
+	    if (dif < 0.0)
+		next
+	    sumi = sumi + dif
+	    sumix = sumix + i * dif
+	}
+
+	if (fp_equalr (sumi, 0.0))
+	    shift = (1.0 + npts) / 2.0
+	else
+	    shift = sumix / sumi
+end
+
+
+define	MIN_WIDTH	3.		# minimum centering width
+define	EPSILON		0.001		# accuracy of centering
+define	EPSILON1	0.005		# tolerance for convergence check
+define	ITERATIONS	100		# maximum number of iterations
+define	MAX_DXCHECK	3		# look back for failed convergence
+define	INTERPTYPE	II_SPLINE3	# image interpolation type
+
+
+# RG_X1DCENTER -- Locate the center of a one dimensional feature.
+# A value of INDEF is returned in the centering fails for any reason.
+# This procedure just sets up the data and adjusts for emission or
+# absorption features.  The actual centering is done by C1D_CENTER.
+
+procedure rg_x1dcenter (x, data, npts, xc, ic, width, type, radius, threshold)
+
+real	x				#I initial guess
+real	data[npts]			#I data points
+int	npts				#I number of data points
+real	xc				#O computed center
+real	ic				#O intensity at computed center
+real	width				#I feature width
+int	type				#I feature type
+real	radius				#I centering radius
+real	threshold			#I minimum range in feature
+
+int	x1, x2, nx
+real	a, b, rad, wid
+pointer	sp, data1
+
+begin
+	# Check starting value.
+	if (IS_INDEF(x) || (x < 1) || (x > npts)) {
+	    xc = INDEF
+	    ic = INDEF
+	    return
+	}
+
+	# Set minimum width and error radius.  The minimum in the error radius
+	# is for defining the data window.  The user error radius is used to
+	# check for an error in the derived center at the end of the centering.
+
+	wid = max (width, MIN_WIDTH)
+	rad = max (2., radius)
+
+	# Determine the pixel value range around the initial center, including
+	# the width and error radius buffer.  Check for a minimum range.
+
+	x1 = max (1., x - wid / 2 - rad - wid)
+	x2 = min (real (npts), x + wid / 2 + rad + wid + 1)
+	nx = x2 - x1 + 1
+	call alimr (data[x1], nx, a, b)
+	if (b - a < threshold) {
+	    xc = INDEF
+	    ic = INDEF
+	    return
+	}
+
+	# Allocate memory for the continuum subtracted data vector.  The X
+	# range is just large enough to include the error radius and the
+	# half width.
+
+	x1 = max (1., x - wid / 2 - rad)
+	x2 = min (real (npts), x + wid / 2 + rad + 1)
+	nx = x2 - x1 + 1
+
+	call smark (sp)
+	call salloc (data1, nx, TY_REAL)
+	call amovr (data[x1], Memr[data1], nx)
+
+	# Make the centering data positive, subtract the continuum, and
+	# apply a threshold to eliminate noise spikes.
+
+	switch (type) {
+	case EMISSION:
+	    a = 0.
+	    call asubkr (data[x1], a + threshold, Memr[data1], nx)
+	    call amaxkr (Memr[data1], 0., Memr[data1], nx)
+	case ABSORPTION:
+	    call anegr (data[x1], Memr[data1], nx)
+	    call asubkr (Memr[data1], threshold - b, Memr[data1], nx)
+	    call amaxkr (Memr[data1], 0., Memr[data1], nx)
+	default:
+	    call error (0, "Unknown feature type")
+	}
+
+	# Determine the center.
+	call rg_xcenter (x - x1 + 1, Memr[data1], nx, xc, ic, wid)
+
+	# Check user centering error radius.
+	if (!IS_INDEF(xc)) {
+	    xc = xc + x1 - 1
+	    if (abs (x - xc) > radius) {
+		xc = INDEF
+		ic = INDEF
+	    }
+	}
+
+	# Free memory and return the center position.
+	call sfree (sp)
+end
+
+
+# RG_XCENTER -- One dimensional centering algorithm.
+
+procedure rg_xcenter (x, data, npts, xc, ic, width)
+
+real	x				#I starting guess
+int	npts				#I number of points in data vector
+real	data[npts]			#I data vector
+real	xc				#O computed xc
+real	ic				#O computed intensity at xc
+real	width				#I centering width
+
+int	i, j, iteration, dxcheck
+real	hwidth, dx, dxabs, dxlast
+real	a, b, sum1, sum2, intgrl1, intgrl2
+pointer	asi1, asi2, sp, data1
+
+real	asigrl(), asieval()
+
+define	done_	99
+
+begin
+	# Find the nearest local maxima as the starting point.
+	# This is required because the threshold limit may have set
+	# large regions of the data to zero and without a gradient
+	# the centering will fail.
+
+	i = x
+	for (i=x+.5; (i<npts) && (data[i]<=data[i+1]); i=i+1)
+	    ;
+	for (j=x+.5; (j>1) && (data[j]<=data[j-1]); j=j-1)
+	    ;
+
+	if (i-x < x-j)
+	    xc = i 
+	else
+	    xc = j
+
+	# Check data range.
+	hwidth = width / 2
+	if ((xc - hwidth < 1) || (xc + hwidth > npts)) {
+	    xc = INDEF
+	    ic = INDEF
+	    return
+	}
+
+	# Set interpolation functions.
+	call asiinit (asi1, INTERPTYPE)
+	call asiinit (asi2, INTERPTYPE)
+	call asifit (asi1, data, npts)
+
+	# Allocate, compute, and interpolate the x*y values.
+	call smark (sp)
+	call salloc (data1, npts, TY_REAL)
+	do i = 1, npts
+	    Memr[data1+i-1] = data[i] * i
+	call asifit (asi2, Memr[data1], npts)
+	call sfree (sp)
+
+	# Iterate to find center.  This loop exits when 1) the maximum
+	# number of iterations is reached, 2) the delta is less than
+	# the required accuracy (criterion for finding a center), 3)
+	# there is a problem in the computation, 4) successive steps
+	# continue to exceed the minimum delta.
+
+	dxlast = 1.
+	do iteration = 1, ITERATIONS {
+
+	    # Triangle centering function.
+	    a = xc - hwidth
+	    b = xc - hwidth / 2
+	    intgrl1 = asigrl (asi1, a, b)
+	    intgrl2 = asigrl (asi2, a, b)
+	    sum1 = (xc - hwidth) * intgrl1 - intgrl2
+	    sum2 = -intgrl1
+	    a = b
+	    b = xc + hwidth / 2
+	    intgrl1 = asigrl (asi1, a, b)
+	    intgrl2 = asigrl (asi2, a, b)
+	    sum1 = sum1 - xc * intgrl1 + intgrl2
+	    sum2 = sum2 + intgrl1
+	    a = b
+	    b = xc + hwidth
+	    intgrl1 = asigrl (asi1, a, b)
+	    intgrl2 = asigrl (asi2, a, b)
+	    sum1 = sum1 + (xc + hwidth) * intgrl1 - intgrl2
+	    sum2 = sum2 - intgrl1
+
+	    # Return no center if sum2 is zero.
+	    if (sum2 == 0.)
+		break
+
+	    # Limit dx change in one iteration to 1 pixel.
+	    dx = max (-1., min (1., sum1 / abs (sum2)))
+	    dxabs = abs (dx)
+	    xc = xc + dx
+	    ic = asieval (asi1, xc)
+
+	    # Check data range.  Return no center if at edge of data.
+	    if ((xc - hwidth < 1) || (xc + hwidth > npts))
+		break
+
+	    # Convergence tests.
+	    if (dxabs < EPSILON)
+		goto done_
+	    if (dxabs > dxlast + EPSILON1) {
+		dxcheck = dxcheck + 1
+		if (dxcheck > MAX_DXCHECK)
+		    break
+	    } else {
+		dxcheck = 0
+	        dxlast = dxabs
+	    }
+	}
+
+	# If we get here then no center was found.
+	xc = INDEF
+	ic = INDEF
+
+done_	call asifree (asi1)
+	call asifree (asi2)
+end
+
+
+# RG_IPARAB -- Compute the coefficients of the parabola through three
+# evenly spaced points.
+
+procedure rg_iparab (x, y, c)
+
+real	x[NPARS_PARABOLA]		#I input x values
+real	y[NPARS_PARABOLA]		#I input y values
+real	c[NPARS_PARABOLA]		#O computed coefficients
+
+begin
+	c[3] = (y[1]-y[2]) * (x[2]-x[3]) / (x[1]-x[2]) - (y[2]-y[3])
+        c[3] = c[3] / ((x[1]**2-x[2]**2) * (x[2]-x[3]) / (x[1]-x[2]) -
+                (x[2]**2-x[3]**2))
+
+        c[2] = (y[1] - y[2]) - c[3] * (x[1]**2 - x[2]**2)
+        c[2] = c[2] / (x[1] - x[2])
+
+        c[1] = y[1] - c[2] * x[1] - c[3] * x[1]**2
+end
+
+
+# RG_POLYFIT -- Evaluate an nth order polynomial.
+
+procedure rg_polyfit (x, nvars, p, np, z)
+
+real	x		#I position coordinate
+int	nvars		#I number of variables
+real	p[ARB]		#I coefficients of polynomial
+int	np		#I number of parameters 
+real	z		#O function return
+
+int	i
+real	r
+
+begin
+	r = 0.0
+	do i = 2, np 
+	    r = r + x**(i-1) * p[i]
+	z = p[1] + r
+end
+
+
+# RG_DPOLYFIT -- Evaluate an nth order polynomial and its derivatives.
+
+procedure rg_dpolyfit (x, nvars, p, dp, np, z, der)
+
+real	x		#I position coordinate
+int	nvars		#I number of variables
+real	p[ARB]		#I coefficients of polynomial
+real	dp[ARB]		#I parameter derivative increments
+int	np		#I number of parameters
+real	z		#O function value
+real	der[ARB]	#O derivatives
+
+int	i
+
+begin
+	der[1] = 1.0
+	z = 0.0
+	do i = 2, np {
+	    der[i] = x ** (i-1)
+	    z = z + x**(i-1) * p[i]  
+	}
+	z = p[1] + z
+end
diff --git a/pkg/images/immatch/src/xregister/rgxgpars.x b/pkg/images/immatch/src/xregister/rgxgpars.x
new file mode 100644
index 00000000..82943730
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxgpars.x
@@ -0,0 +1,68 @@
+include "xregister.h"
+
+# RG_XGPARS -- Read in the XREGISTER task algorithm parameters.
+
+procedure rg_xgpars (xc)
+
+pointer	xc		#I pointer to the main structure
+
+int	xlag, ylag, xwindow, ywindow, xcbox, ycbox
+pointer	sp, str
+int	clgwrd(), clgeti()
+real	clgetr()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Initialize the correlation structure.
+	call rg_xinit (xc, clgwrd ("correlation", Memc[str], SZ_LINE,
+	    XC_CTYPES))
+
+	# Fetch the initial shift information.
+	xlag = clgeti ("xlag")
+	ylag = clgeti ("ylag")
+	call rg_xseti (xc, IXLAG, xlag)
+	call rg_xseti (xc, IYLAG, ylag)
+	call rg_xseti (xc, XLAG, xlag)
+	call rg_xseti (xc, YLAG, ylag)
+	call rg_xseti (xc, DXLAG, clgeti ("dxlag"))
+	call rg_xseti (xc, DYLAG, clgeti ("dylag"))
+
+	# Get the background value computation parameters.
+	call rg_xseti (xc, BACKGRD, clgwrd ("background", Memc[str], SZ_LINE,
+	    XC_BTYPES))
+	call rg_xsets (xc, BSTRING, Memc[str])
+	call rg_xseti (xc, BORDER, clgeti ("border"))
+	call rg_xsetr (xc, LOREJECT, clgetr ("loreject"))
+	call rg_xsetr (xc, HIREJECT, clgetr ("hireject"))
+	call rg_xsetr (xc, APODIZE, clgetr ("apodize"))
+	call rg_xseti (xc, FILTER, clgwrd ("filter", Memc[str], SZ_LINE,
+	    XC_FTYPES))
+	call rg_xsets (xc, FSTRING, Memc[str])
+
+	# Get the window parameters and force the window size to be odd.
+	xwindow = clgeti ("xwindow")
+	if (mod (xwindow,2) == 0)
+	    xwindow = xwindow + 1
+	call rg_xseti (xc, XWINDOW, xwindow)
+	ywindow = clgeti ("ywindow")
+	if (mod (ywindow,2) == 0)
+	    ywindow = ywindow + 1
+	call rg_xseti (xc, YWINDOW, ywindow)
+
+	# Get the peak fitting parameters.
+	call rg_xseti (xc, PFUNC, clgwrd ("function", Memc[str], SZ_LINE,
+	    XC_PTYPES))
+	xcbox = clgeti ("xcbox")
+	if (mod (xcbox,2) == 0)
+	    xcbox = xcbox + 1
+	call rg_xseti (xc, XCBOX, xcbox)
+	ycbox = clgeti ("ycbox")
+	if (mod (ycbox,2) == 0)
+	    ycbox = ycbox + 1
+	call rg_xseti (xc, YCBOX, ycbox)
+
+	call sfree (sp)
+end
diff --git a/pkg/images/immatch/src/xregister/rgxicorr.x b/pkg/images/immatch/src/xregister/rgxicorr.x
new file mode 100644
index 00000000..e96c6dec
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxicorr.x
@@ -0,0 +1,583 @@
+include <imhdr.h>
+include <fset.h>
+include <ctype.h>
+include "xregister.h"
+
+define	HELPFILE	"immatch$src/xregister/xregister.key"
+define	OHELPFILE	"immatch$src/xregister/oxregister.key"
+
+define	XC_PCONTOUR	1
+define	XC_PLINE	2
+define	XC_PCOL		3
+
+
+# RG_XICORR -- Compute the shifts for each image interactively using
+# cross-correlation techniques.
+
+int procedure rg_xicorr (imr, im1, im2, db, dformat, reglist, tfd, xc, gd, id)
+
+pointer	imr		#I/O pointer to the reference image
+pointer	im1		#I/O pointer to the input image
+pointer	im2		#I/O pointer to the output image
+pointer	db		#I/O pointer to the shifts database file
+int	dformat		#I is the shifts file in database format
+int	reglist		#I/O the regions list descriptor
+int	tfd		#I/O the transform file descriptor
+pointer	xc		#I pointer to the cross-corrrelation structure
+pointer	gd		#I the graphics stream pointer
+pointer	id		#I the display stream pointer
+
+int	newdata, newcross, newcenter, wcs, key, cplottype, newplot
+int	ip, ncolr, nliner
+pointer	sp, cmd
+real	xshift, yshift, wx, wy
+int	rg_xstati(), rg_icross(), clgcur(), rg_xgtverify(), rg_xgqverify()
+int	ctoi()
+pointer	rg_xstatp()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	# Initialize.
+	newdata = YES
+	newcross = YES
+	newcenter = YES
+	ncolr = (1 + rg_xstati (xc, XWINDOW)) / 2
+	nliner = (1 + rg_xstati (xc, YWINDOW)) / 2
+	cplottype = XC_PCONTOUR 
+	newplot = YES
+	xshift = 0.0
+	yshift = 0.0
+
+	# Compute the cross-correlation function for the first region
+	# and print the results.
+	if (rg_xstati (xc, NREGIONS) <= 0) {
+	    call gclear (gd)
+	    call printf ("The regions list is empty\n")
+	} else if (rg_icross (xc, imr, im1, rg_xstati (xc, CREGION)) != ERR) {
+	    call rg_xcplot (xc, gd, ncolr, nliner, cplottype)
+	    call rg_fit (xc, rg_xstati (xc, CREGION), gd, xshift, yshift)
+	    call rg_xpwrec (xc, rg_xstati (xc, CREGION))
+	    newdata = NO
+	    newcross = NO
+	    newcenter = NO
+	    newplot = NO
+	} else {
+	    call gclear (gd)
+	    call printf (
+	        "Error computing X-correlation function for region %d\n")
+		call pargi (rg_xstati (xc, CREGION))
+	}
+
+
+	# Loop over the cursor commands.
+	while (clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) !=
+	    EOF) {
+
+	    switch (key) {
+
+	    # Print the help page.
+	    case '?':
+	        call gpagefile (gd, HELPFILE, "")
+
+	    # Redraw the current plot.
+	    case 'r':
+		newplot = YES
+
+	    # Draw a contour plot of the cross-correlation function.
+	    case 'c':
+		if (cplottype != XC_PCONTOUR)
+		    newplot = YES
+		ncolr = (rg_xstati (xc, XWINDOW) + 1) / 2
+		nliner = (rg_xstati (xc, YWINDOW) + 1) / 2
+		cplottype = XC_PCONTOUR 
+
+	    # Plot a column of the cross-correlation function.
+	    case 'x':
+		if (cplottype != XC_PCOL)
+		    newplot = YES
+		if (cplottype == XC_PCONTOUR) {
+		    ncolr = nint (wx)
+		    nliner = nint (wy)
+		} else if (cplottype == XC_PLINE) {
+		    ncolr = nint (wx)
+		} 
+		cplottype = XC_PCOL 
+
+	    # Plot a line of the cross-correlation function.
+	    case 'y':
+		if (cplottype != XC_PLINE)
+		    newplot = YES
+		if (cplottype == XC_PCONTOUR) {
+		    ncolr = nint (wx)
+		    nliner = nint (wy)
+		} else if (cplottype == XC_PCOL) {
+		    ncolr = nint (wx)
+		} 
+		cplottype = XC_PLINE 
+
+	    # Quit the task gracefully.
+	    case 'q':
+		if (rg_xgqverify ("xregister", db, dformat, xc, key) == YES) {
+		    call sfree (sp)
+		    return (rg_xgtverify (key))
+		}
+
+	    # The Data overlay menu.
+	    case 'o':
+		#call gdeactivate (gd, 0)
+		call rg_xoverlay (gd, xc, rg_xstati (xc, CREGION), imr, im1)
+		#call greactivate (gd, 0)
+		newplot = YES
+
+	    # Process colon commands.
+	    case ':':
+		for (ip = 1; IS_WHITE(Memc[cmd+ip-1]); ip = ip + 1)
+		    ;
+		switch (Memc[cmd+ip-1]) {
+		case 'x':
+		    if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') {
+	                call rg_xcolon (gd, xc, imr, im1, im2, db, dformat,
+			    tfd, reglist, Memc[cmd], newdata, newcross,
+			    newcenter)
+		    } else {
+			ip = ip + 1
+			if (ctoi (Memc[cmd], ip, ncolr) <= 0)
+			    ncolr = (1 + rg_xstati (xc, XWINDOW)) / 2
+			cplottype = XC_PCOL
+			newplot = YES
+		    }
+		case 'y':
+		    if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') {
+	                call rg_xcolon (gd, xc, imr, im1, im2, db, dformat,
+			    tfd, reglist, Memc[cmd], newdata, newcross,
+			    newcenter)
+		    } else {
+			ip = ip + 1
+			if (ctoi (Memc[cmd], ip, nliner) <= 0)
+			    nliner = (1 + rg_xstati (xc, YWINDOW)) / 2
+			cplottype = XC_PLINE
+			newplot = YES
+		    }
+		default:
+	            call rg_xcolon (gd, xc, imr, im1, im2, db, dformat, tfd,
+		        reglist, Memc[cmd], newdata, newcross, newcenter)
+		}
+
+	    # Compute an image lag interactively.
+	    case 't':
+		call gdeactivate (gd, 0)
+		call rg_itransform (xc, imr, im1, id)
+		newdata = YES; newcross = YES; newcenter = YES
+		call greactivate (gd, 0)
+
+	     # Write the parameters to the parameter file.
+	     case 'w':
+		call rg_pxpars (xc)
+
+	    case 'f':
+
+		if (rg_xstati (xc, NREGIONS) > 0) {
+
+	    	    if (newdata == YES) {
+			call rg_xcindefr (xc, rg_xstati(xc,CREGION))
+			newdata = NO
+	    	    }
+
+	    	    if (newcross == YES) {
+			call printf (
+			    "Recomputing X-correlation function ...\n")
+			if (rg_icross (xc, imr, im1, rg_xstati (xc,
+			        CREGION)) != ERR) {
+			    ncolr = (1 + rg_xstati (xc, XWINDOW)) / 2
+			    if (IM_NDIM(imr) == 1)
+				nliner = 1
+			    else
+			        nliner = (1 + rg_xstati (xc, YWINDOW)) / 2
+	    	            call rg_xcplot (xc, gd, ncolr, nliner, cplottype)
+	    	            call rg_fit (xc, rg_xstati (xc, CREGION), gd,
+				xshift, yshift)
+	    	            call rg_xpwrec (xc, rg_xstati (xc, CREGION))
+	    	            newcross = NO
+	    	            newcenter = NO
+	    	            newplot = NO
+		        } else {
+	    	            call printf (
+	               "Error computing X-correlation function for region %d\n")
+			        call pargi (rg_xstati (xc, CREGION))
+		        }
+	            }
+
+	    	    if (newcenter == YES) {
+	    		call rg_fit (xc, rg_xstati (xc, CREGION), gd,
+			    xshift, yshift)
+	    		call rg_xpwrec (xc, rg_xstati (xc, CREGION))
+	    		newcenter = NO
+	    	    }
+
+		} else
+		    call printf ("The regions list is empty\n")
+
+
+
+	    # Do nothing gracefully.
+	    default:
+		call printf ("Unknown or ambiguous keystroke command\n")
+	    }
+
+	    # Replot the correlation function.
+	    if (newplot == YES) {
+	        if (newdata == YES) {
+		    call printf (
+		        "Warning: X-correlation function should be refit\n")
+	        } else if (newcross == YES) {
+		    call printf (
+		        "Warning: X-correlation function should be refit\n")
+	        } else if (newcenter == YES) {
+		    call printf (
+		        "Warning: X-correlation function should be refit\n")
+	        } else if (rg_xstatp (xc, XCOR) != NULL) {
+	            call rg_xcplot (xc, gd, ncolr, nliner, cplottype)
+	            call rg_xpwrec (xc, rg_xstati (xc, CREGION))
+		    newplot = NO
+	        } else {
+		    call printf (
+		        "Warning: X-correlation function is undefined\n")
+	        }
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# RG_XOVERLAY -- The image overlay plot menu.
+
+procedure rg_xoverlay (gd, xc, nreg, imr, im1)
+
+pointer	gd		#I graphics stream pointer
+pointer	xc		#I pointer to the crosscor structure
+int	nreg		#I the current region number
+pointer	imr		#I pointer to the reference image
+pointer	im1		#I pointer to the input image
+
+int	ip, wcs, key, ixlag, iylag, ixshift, iyshift
+int	nrimcols, nrimlines, nimcols, nimlines, ncolr, ncoli, nliner, nlinei
+pointer	sp, cmd
+real	wx, wy, rxlag, rylag, xshift, yshift
+int	clgcur(), ctoi(), rg_xstati()
+pointer	rg_xstatp()
+
+begin
+	if (gd == NULL)
+	    return
+
+	nrimcols = IM_LEN(imr,1)
+	if (IM_NDIM(imr) == 1)
+	    nrimlines = 1
+	else
+	    nrimlines = IM_LEN(imr,2)
+	nimcols = IM_LEN(im1,1)
+	if (IM_NDIM(im1) == 1)
+	    nimlines = 1
+	else
+	    nimlines = IM_LEN(im1,2)
+	if (rg_xstati (xc, NREFPTS) > 0) {
+	    wx = (1. + nrimcols) / 2.0
+	    wy = (1. + nrimlines) / 2.0
+	    call rg_etransform (xc, wx, wy, rxlag, rylag)
+	    ixlag = rxlag - wx
+	    iylag = rylag - wy
+	} else {
+	    ixlag = rg_xstati (xc, XLAG) 
+	    iylag = rg_xstati (xc, YLAG) 
+	}
+	xshift = -Memr[rg_xstatp(xc,XSHIFTS)+nreg-1]
+	yshift = -Memr[rg_xstatp(xc,YSHIFTS)+nreg-1]
+
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	while (clgcur ("icommands", wx, wy, wcs, key, Memc[cmd],
+	    SZ_LINE) != EOF) {
+
+	    switch (key) {
+
+	    # Print the help menu.
+	    case '?':
+		call gdeactivate (gd, 0)
+		call pagefile (OHELPFILE, "")
+		call greactivate (gd, 0)
+
+	    # Quit.
+	    case 'q':
+		break
+
+	    # Plot the same line of the reference and input image.
+	    case 'l':
+		call rg_xpline (gd, imr, im1, nint (wy), 0, 0)
+
+	    # Plot the same column of the reference and input image
+	    case 'c':
+		call rg_xpcol (gd, imr, im1, nint (wx), 0, 0)
+
+	    case 'y':
+		call rg_xpline (gd, imr, im1, nint (wy), ixlag, iylag)
+
+	    case 'x':
+		call rg_xpcol (gd, imr, im1, nint (wx), ixlag, iylag)
+
+	    case 'h':
+		call rg_xpline (gd, imr, im1, nint (wy), nint (xshift),
+		    nint (yshift))
+
+	    case 'v':
+		call rg_xpcol (gd, imr, im1, nint (wx), nint (xshift),
+		    nint (yshift))
+
+	    case ':':
+		ip = 1
+	        call rg_cokeys (Memc[cmd], ip, SZ_LINE, key)
+		switch (key) {
+	        case 'l':
+		    ixshift = 0
+		    if (ctoi (Memc[cmd], ip, nliner) <= 0)
+			nliner = (1 + nrimlines) / 2 
+		    nliner = max (1, min (nliner, nrimlines))
+		    if (ctoi (Memc[cmd], ip, nlinei) <= 0)
+			nlinei = nliner
+		    iyshift = nlinei - nliner
+		    call rg_xpline (gd, imr, im1, nliner, ixshift, iyshift)
+
+		case 'c':
+		    if (ctoi (Memc[cmd], ip, ncolr) <= 0)
+		        ncolr = (1 + nrimcols) / 2
+		    ncolr = max (1, min (ncolr, nrimcols))
+		    if (ctoi (Memc[cmd], ip, ncoli) <= 0)
+			ncoli = ncolr
+		    ncoli = max (1, min (ncoli, nimcols))
+		    ixshift = ncoli - ncolr
+		    iyshift = 0
+		    call rg_xpcol (gd, imr, im1, ncolr, ixshift, iyshift)
+
+		case 'y':
+		    if (ctoi (Memc[cmd], ip, nliner) <= 0)
+		 	nliner = (1 + nrimlines) / 2
+		    nliner = max (1, min (nliner, nrimlines))
+		    call rg_xpline (gd, imr, im1, nliner, ixlag, iylag)
+
+		case 'x':
+		    if (ctoi (Memc[cmd], ip, ncolr) <= 0)
+			ncolr = (1 + nrimcols) / 2
+		    ncolr = max (1, min (ncolr, nrimcols))
+		    call rg_xpcol (gd, imr, im1, ncolr, ixlag, iylag)
+
+		case 'h':
+		    if (ctoi (Memc[cmd], ip, nliner) <= 0)
+		 	nliner = (1 + nrimlines) / 2
+		    nliner = max (1, min (nliner, nrimlines))
+		    call rg_xpline (gd, imr, im1, nliner, nint (xshift),
+		        nint (yshift))
+
+		case 'v':
+		    if (ctoi (Memc[cmd], ip, ncolr) <= 0)
+			ncolr = (1 + nrimcols) / 2
+		    ncolr = max (1, min (ncolr, nrimcols))
+		    call rg_xpcol (gd, imr, im1, ncolr, nint (xshift),
+		        nint (yshift))
+		default:
+		    call printf ("Ambiguous or unknown overlay menu command\n")
+		}
+	    case 'g':
+		while (clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd],
+		    SZ_LINE) != EOF) {
+		    if (key == 'q')
+			break
+		}
+	    default:
+		call printf ("Ambiguous or unknown overlay menu command\n")
+	    }
+
+	}
+
+	call sfree (sp)
+end
+
+
+# RG_XCPLOT -- Draw the default plot of the cross-correlation function.
+
+procedure rg_xcplot (xc, gd, col, line, plottype)
+
+pointer	xc		#I pointer to cross-correlation structure
+pointer	gd		#I pointer to the graphics stream
+int	col		#I column of cross-correlation function to plot
+int	line		#I line of cross-correlation function to plot
+int	plottype	#I the default plot type
+
+int	nreg, xwindow, ywindow
+pointer	sp, title, str, prc1, prc2, prl1, prl2
+int	rg_xstati(), strlen()
+pointer	rg_xstatp()
+
+begin
+	if (gd == NULL)
+	    return
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (title, SZ_LINE, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Get the regions.
+	nreg = rg_xstati (xc, CREGION)
+	prc1 = rg_xstatp (xc, RC1)
+	prc2 = rg_xstatp (xc, RC2)
+	prl1 = rg_xstatp (xc, RL1)
+	prl2 = rg_xstatp (xc, RL2)
+
+	# Initialize the window size.
+	xwindow = rg_xstati (xc, XWINDOW)
+	if ((Memi[prl2+nreg-1] - Memi[prl1+nreg-1] + 1) == 1)
+	    ywindow = 1
+	else
+	    ywindow = rg_xstati (xc, YWINDOW)
+
+	# Construct a title.
+	call sprintf (Memc[title], SZ_LINE,
+	    "Reference: %s  Image: %s  Region: [%d:%d,%d:%d]")
+	    call rg_xstats (xc, REFIMAGE, Memc[str], SZ_FNAME)
+	    call pargstr (Memc[str])
+	    call rg_xstats (xc, IMAGE, Memc[str], SZ_FNAME)
+	    call pargstr (Memc[str])
+	    call pargi (Memi[prc1+nreg-1])
+	    call pargi (Memi[prc2+nreg-1])
+	    call pargi (Memi[prl1+nreg-1])
+	    call pargi (Memi[prl2+nreg-1])
+
+	# Draw the plot.
+	if (ywindow == 1) {
+	    call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+	        "\nX-Correlation Function: line %d")
+		call pargi (1)
+	    call rg_xcpline (gd, Memc[title], Memr[rg_xstatp(xc,XCOR)],
+	        xwindow, ywindow, 1)
+	} else {
+	    switch (plottype) {
+	    case XC_PCONTOUR:
+	        call rg_contour (gd, "X-Correlation Function", Memc[title],
+		    Memr[rg_xstatp (xc, XCOR)], xwindow, ywindow)
+	    case XC_PLINE:
+	        call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+	            "\nX-Correlation Function: line %d")
+		    call pargi (line)
+	        call rg_xcpline (gd, Memc[title], Memr[rg_xstatp(xc,XCOR)],
+	            xwindow, ywindow, line)
+	    case XC_PCOL:
+	        call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+	            "\nX-Correlation Function: column %d")
+		    call pargi (col)
+	        call rg_xcpcol (gd, Memc[title], Memr[rg_xstatp(xc,XCOR)],
+	            xwindow, ywindow, col)
+	    default:
+	        call rg_contour (gd, "X-Correlation Function", Memc[title],
+		    Memr[rg_xstatp (xc, XCOR)], xwindow, ywindow)
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# RG_COKEYS -- Fetch the first keystroke of a colon command.
+
+procedure rg_cokeys (cmd, ip, maxch, key)
+
+char	cmd[ARB]		#I the command string
+int	ip			#I/O pointer into the command string
+int	maxch			#I maximum number of characters
+int	key			#O the keystroke
+
+begin
+	ip = 1
+	while (IS_WHITE(cmd[ip]) && cmd[ip] != EOS && ip <= maxch)
+	    ip = ip + 1
+
+	if (cmd[ip] == EOS && ip > maxch)
+	    key = EOS
+	else {
+	    key = cmd[ip]
+	    ip = ip + 1
+	}
+end
+
+
+define	QUERY "Hit [return=continue, n=next image, q=quit, w=quit and update parameters]: "
+
+# RG_XGQVERIFY -- Print a message on the status line asking the user if they
+# really want to quit, returning YES if they really want to quit, NO otherwise.
+
+int procedure rg_xgqverify (task, db, dformat, rg, ch)
+
+char	task[ARB]	#I the calling task name
+pointer	db		#I pointer to the shifts database file
+int	dformat		#I is the shifts file in database format
+pointer	rg		#I pointer to the task structure
+int	ch		#I the input keystroke command
+
+int	wcs, stat
+pointer	sp, cmd
+real	wx, wy
+bool	streq()
+int	clgcur()
+
+begin
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	# Print the status line query in reverse video and get the keystroke.
+	call printf (QUERY)
+	#call flush (STDOUT)
+	if (clgcur ("gcommands", wx, wy, wcs, ch, Memc[cmd], SZ_LINE) == EOF)
+	    ;
+
+	# Process the command.
+	if (ch == 'q') {
+	    call rg_xwrec (db, dformat, rg)
+	    stat = YES
+	} else if (ch == 'w') {
+	    call rg_xwrec (db, dformat, rg)
+	    if (streq ("xregister", task))
+		call rg_pxpars (rg)
+	    stat = YES
+	} else if (ch == 'n') {
+	    call rg_xwrec (db, dformat, rg)
+	    stat = YES
+	} else {
+	    stat = NO
+	}
+
+	call sfree (sp)
+	return (stat)
+end
+
+
+# RG_XGTVERIFY -- Verify whether or not the user truly wishes to quit the
+# task.
+
+int procedure rg_xgtverify (ch)
+
+int	ch		#I the input keystroke command
+
+begin
+	if (ch == 'q') {
+	    return (YES)
+	} else if (ch == 'w') {
+	    return (YES)
+	} else if (ch == 'n') {
+	    return (NO)
+	} else {
+	    return (NO)
+	}
+end
diff --git a/pkg/images/immatch/src/xregister/rgximshift.x b/pkg/images/immatch/src/xregister/rgximshift.x
new file mode 100644
index 00000000..08cb3f62
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgximshift.x
@@ -0,0 +1,391 @@
+include <imhdr.h>
+include <imset.h>
+include <math/iminterp.h>
+
+define	NYOUT		16	# number of lines output at once
+define	NMARGIN		3	# number of boundary pixels required	
+define	NMARGIN_SPLINE3	16	# number of spline boundary pixels required
+
+
+# RG_XSHIFTIM - Shift a 1 or 2D image by a fractional pixel amount
+# x and y
+
+procedure rg_xshiftim (im1, im2, xshift, yshift, interpstr, boundary_type,
+        constant)
+
+pointer		im1		#I pointer to input image
+pointer		im2		#I pointer to output image
+real		xshift		#I shift in x direction
+real		yshift		#I shift in y direction
+char		interpstr[ARB]	#I type of interpolant
+int		boundary_type	#I type of boundary extension
+real		constant	#I value of constant for boundary extension
+
+int	interp_type
+pointer	sp, str
+bool	fp_equalr()
+int	strdic()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+	interp_type = strdic (interpstr, Memc[str], SZ_FNAME, II_BFUNCTIONS)
+
+	if (interp_type == II_NEAREST)
+	    call rg_xishiftim (im1, im2, nint (xshift), nint (yshift),
+	        interp_type, boundary_type, constant)
+	else if (fp_equalr (xshift, real (int (xshift))) && fp_equalr (yshift,
+	    real (int (xshift))))
+	    call rg_xishiftim (im1, im2, int (xshift), int (yshift),
+	        interp_type, boundary_type, constant)
+	else
+	    call rg_xfshiftim (im1, im2, xshift, yshift, interpstr,
+		boundary_type, constant)
+	call sfree (sp)
+end
+
+
+# RG_XISHIFTIM -- Shift a 2-D image by integral pixels in x and y.
+
+procedure rg_xishiftim (im1, im2, nxshift, nyshift, interp_type, boundary_type,
+	constant)
+
+pointer	im1			#I pointer to the input image
+pointer	im2			#I pointer to the output image
+int	nxshift, nyshift	#I shift in x and y
+int	interp_type		#I type of interpolant
+int	boundary_type		#I type of boundary extension
+real	constant		#I constant for boundary extension
+
+int	ixshift, iyshift
+pointer	buf1, buf2
+long	v[IM_MAXDIM]
+int	ncols, nlines, nbpix
+int	i, x1col, x2col, yline
+
+int	impnls(), impnli(), impnll(), impnlr(), impnld(), impnlx()
+pointer	imgs2s(), imgs2i(), imgs2l(), imgs2r(), imgs2d(), imgs2x()
+errchk	impnls, impnli, impnll, impnlr, impnld, impnlx
+errchk	imgs2s, imgs2i, imgs2l, imgs2r, imgs2d, imgs2x
+string	wrerr "ISHIFTXY: Error writing in image."
+
+begin
+	ixshift = nxshift
+	iyshift = nyshift
+
+	ncols = IM_LEN(im1,1)
+	nlines = IM_LEN(im1,2)
+
+	# Cannot shift off image.
+	if (ixshift < -ncols || ixshift > ncols)
+	    call error (3, "ISHIFTXY: X shift out of bounds.")
+	if (iyshift < -nlines || iyshift > nlines)
+	    call error (4, "ISHIFTXY: Y shift out of bounds.")
+
+	# Calculate the shift.
+	switch (boundary_type) {
+	case BT_CONSTANT,BT_REFLECT,BT_NEAREST:
+	    ixshift = min (ncols, max (-ncols, ixshift))
+	    iyshift = min (nlines, max (-nlines, iyshift))
+	case BT_WRAP:
+	    ixshift = mod (ixshift, ncols)
+	    iyshift = mod (iyshift, nlines)
+	}
+
+	# Set the boundary extension values.
+	nbpix = max (abs (ixshift), abs (iyshift))
+	call imseti (im1, IM_NBNDRYPIX, nbpix)
+	call imseti (im1, IM_TYBNDRY, boundary_type)
+	if (boundary_type == BT_CONSTANT)
+	    call imsetr (im1, IM_BNDRYPIXVAL, constant)
+
+	# Get column boundaries in the input image.
+	x1col = max (-ncols + 1, - ixshift + 1) 
+	x2col = min (2 * ncols,  ncols - ixshift)
+
+	call amovkl (long (1), v, IM_MAXDIM)
+
+	# Shift the image using the appropriate data type operators.
+	switch (IM_PIXTYPE(im1)) {
+	case TY_SHORT:
+	    do i = 1, nlines {
+	        if (impnls (im2, buf2, v) == EOF)
+		    call error (5, wrerr)
+		yline = i - iyshift
+		buf1 = imgs2s (im1, x1col, x2col, yline, yline)
+		if (buf1 == EOF)
+		    call error (5, wrerr)
+		call amovs (Mems[buf1], Mems[buf2], ncols)
+	    }
+	case TY_INT:
+	    do i = 1, nlines {
+	        if (impnli (im2, buf2, v) == EOF)
+		    call error (5, wrerr)
+		yline = i - iyshift
+		buf1 = imgs2i (im1, x1col, x2col, yline, yline)
+		if (buf1 == EOF)
+		    call error (5, wrerr)
+		call amovi (Memi[buf1], Memi[buf2], ncols)
+	    }
+	case TY_USHORT, TY_LONG:
+	    do i = 1, nlines {
+	        if (impnll (im2, buf2, v) == EOF)
+		    call error (5, wrerr)
+		yline = i - iyshift
+		buf1 = imgs2l (im1, x1col, x2col, yline, yline)
+		if (buf1 == EOF)
+		    call error (5, wrerr)
+		call amovl (Meml[buf1], Meml[buf2], ncols)
+	    }
+	case TY_REAL:
+	    do i = 1, nlines {
+	        if (impnlr (im2, buf2, v) == EOF)
+		    call error (5, wrerr)
+		yline = i - iyshift
+		buf1 = imgs2r (im1, x1col, x2col, yline, yline)
+		if (buf1 == EOF)
+		    call error (5, wrerr)
+		call amovr (Memr[buf1], Memr[buf2], ncols)
+	    }
+	case TY_DOUBLE:
+	    do i = 1, nlines {
+	        if (impnld (im2, buf2, v) == EOF)
+		    call error (0, wrerr)
+		yline = i - iyshift
+		buf1 = imgs2d (im1, x1col, x2col, yline, yline)
+		if (buf1 == EOF)
+		    call error (0, wrerr)
+		call amovd (Memd[buf1], Memd[buf2], ncols)
+	    }
+	case TY_COMPLEX:
+	    do i = 1, nlines {
+	        if (impnlx (im2, buf2, v) == EOF)
+		    call error (0, wrerr)
+		yline = i - iyshift
+		buf1 = imgs2x (im1, x1col, x2col, yline, yline)
+		if (buf1 == EOF)
+		    call error (0, wrerr)
+		call amovx (Memx[buf1], Memx[buf2], ncols)
+	    }
+	default:
+	    call error (6, "ISHIFTXY: Unknown IRAF type.")
+	}
+end
+
+
+
+# RG_XFSHIFTIM -- Shift a 1 or 2D image by a fractional pixel amount
+# in x and y.
+
+procedure rg_xfshiftim (im1, im2, xshift, yshift, interpstr, boundary_type,
+        constant)
+
+pointer		im1		#I pointer to input image
+pointer		im2		#I pointer to output image
+real		xshift		#I shift in x direction
+real		yshift		#I shift in y direction
+char		interpstr[ARB]	#I type of interpolant
+int		boundary_type	#I type of boundary extension
+real		constant	#I value of constant for boundary extension
+
+int	i, interp_type, nsinc, nincr
+int	ncols, nlines, nbpix, fstline, lstline, nxymargin
+int	cin1, cin2, nxin, lin1, lin2, nyin
+int	lout1, lout2, nyout
+real	xshft, yshft, deltax, deltay, dx, dy, cx, ly
+pointer	sp, x, y, msi, sinbuf, soutbuf
+bool	fp_equalr()
+int	msigeti()
+pointer	imps2r()
+
+errchk	imgs2r, imps2r
+errchk	msiinit, msifree, msifit, msigrid
+errchk	smark, salloc, sfree
+
+begin
+	ncols = IM_LEN(im1,1)
+	nlines = IM_LEN(im1,2)
+
+	# Check for out of bounds shift.
+	if (xshift < -ncols || xshift > ncols)
+	    call error (0, "XC_SHIFTIM: X shift out of bounds.")
+	if (yshift < -nlines || yshift > nlines)
+	    call error (0, "XC_SHIFTIM: Y shift out of bounds.")
+
+	# Get the real shift.
+	if (boundary_type == BT_WRAP) {
+	    xshft = mod (xshift, real (ncols))
+	    yshft = mod (yshift, real (nlines))
+	} else {
+	    xshft = xshift
+	    yshft = yshift
+	}
+
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (x, 2 * ncols, TY_REAL)
+	call salloc (y, 2 * nlines, TY_REAL)
+	sinbuf = NULL
+
+	# Define the x and y interpolation coordinates.
+	dx = abs (xshft - int (xshft))
+	if (fp_equalr (dx, 0.0))
+	    deltax = 0.0
+	else if (xshft > 0.)
+	    deltax = 1. - dx
+	else
+	    deltax = dx
+	dy = abs (yshft - int (yshft))
+	if (fp_equalr (dy, 0.0))
+	    deltay = 0.0
+	else if (yshft > 0.)
+	    deltay = 1. - dy
+	else
+	    deltay = dy
+
+	# Initialize the 2-D interpolation routines.
+	call msitype (interpstr, interp_type, nsinc, nincr, cx)
+	if (interp_type == II_BILSINC || interp_type == II_BISINC)
+	    call msisinit (msi, interp_type, nsinc, 1, 1,
+	        deltax - nint (deltax), deltay - nint (deltay), 0.0)
+	else
+	    call msisinit (msi, interp_type, nsinc, 1, 1, cx, cx, 0.0)
+
+	# Set boundary extension parameters.
+	if (interp_type == II_BISPLINE3)
+	    nxymargin = NMARGIN_SPLINE3
+	else if (interp_type == II_BISINC || interp_type == II_BILSINC)
+	    nxymargin = msigeti (msi, II_MSINSINC)
+	else
+	    nxymargin = NMARGIN
+	nbpix = max (int (abs(xshft)+1.0), int (abs(yshft)+1.0)) + nxymargin
+	call imseti (im1, IM_NBNDRYPIX, nbpix)
+	call imseti (im1, IM_TYBNDRY, boundary_type)
+	if (boundary_type == BT_CONSTANT)
+	    call imsetr (im1, IM_BNDRYPIXVAL, constant)
+
+	# Define the x interpolation coordinates.
+	deltax = deltax + nxymargin
+        if (interp_type == II_BIDRIZZLE) {
+            do i = 1, ncols {
+                Memr[x+2*i-2] = i + deltax - 0.5
+                Memr[x+2*i-1] = i + deltax + 0.5
+            }
+        } else {
+	    do i = 1, ncols
+	        Memr[x+i-1] = i + deltax
+	}
+
+	# Define the y interpolation coordinates.
+	deltay = deltay + nxymargin
+        if (interp_type == II_BIDRIZZLE) {
+            do i = 1, NYOUT {
+                Memr[y+2*i-2] = i + deltay - 0.5
+                Memr[y+2*i-1] = i + deltay + 0.5
+            }
+        } else {
+	    do i = 1, NYOUT
+	        Memr[y+i-1] = i + deltay
+	}
+
+	# Define column range in the input image.
+	cx = 1. - nxymargin - xshft
+	if ((cx <= 0.) &&  (! fp_equalr (dx, 0.0)))
+	    cin1 = int (cx) - 1
+	else
+	    cin1 = int (cx)
+	cin2 = ncols - xshft + nxymargin + 1
+	nxin = cin2 - cin1 + 1
+
+	# Loop over output sections.
+	for (lout1 = 1; lout1 <= nlines; lout1 = lout1 + NYOUT) { 
+
+	    # Define range of output lines.
+	    lout2 = min (lout1 + NYOUT - 1, nlines)
+	    nyout = lout2 - lout1 + 1
+
+	    # Define correspoding range of input lines.
+	    ly = lout1 - nxymargin - yshft
+	    if ((ly <= 0) && (! fp_equalr (dy, 0.0)))
+	        lin1 = int (ly) - 1
+	    else
+		lin1 = int (ly)
+	    lin2 = lout2 - yshft + nxymargin + 1
+	    nyin = lin2 - lin1 + 1
+
+	    # Get appropriate input image section and compute the coefficients.
+	    if ((sinbuf == NULL) || (lin1 < fstline) || (lin2 > lstline)) {
+		fstline = lin1
+		lstline = lin2
+		call rg_buf (im1, cin1, cin2, lin1, lin2, sinbuf)
+	        call msifit (msi, Memr[sinbuf], nxin, nyin, nxin)
+	    }
+
+	    # Output the image section.
+	    soutbuf = imps2r (im2, 1, ncols, lout1, lout2)
+	    if (soutbuf == EOF)
+		call error (0, "GSHIFTXY: Error writing output image.")
+
+	    # Evaluate the interpolant.
+	    call msigrid (msi, Memr[x], Memr[y], Memr[soutbuf], ncols, nyout,
+		ncols)
+	}
+
+	call msifree (msi)
+	call sfree (sp)
+end
+
+
+# RG_BUF -- Procedure to provide a buffer of image lines with minimum reads
+
+procedure rg_buf (im, col1, col2, line1, line2, buf)
+
+pointer	im		#I pointer to input image
+int	col1, col2	#I column range of input buffer
+int	line1, line2	#I line range of input buffer
+pointer	buf		#I buffer
+
+int	i, ncols, nlines, nclast, llast1, llast2, nllast
+pointer	buf1, buf2
+
+pointer	imgs2r()
+
+begin
+	ncols = col2 - col1 + 1
+	nlines = line2 - line1 + 1
+
+	if (buf == NULL) {
+	    call malloc (buf, ncols * nlines, TY_REAL)
+	    llast1 = line1 - nlines
+	    llast2 = line2 - nlines
+	} else if ((nlines != nllast) || (ncols != nclast)) {
+	    call realloc (buf, ncols * nlines, TY_REAL)
+	    llast1 = line1 - nlines
+	    llast2 = line2 - nlines
+	}
+
+	if (line1 < llast1) {
+	    do i = line2, line1, -1 {
+		if (i > llast1)
+		    buf1 = buf + (i - llast1) * ncols
+		else
+		    buf1 = imgs2r (im, col1, col2, i, i)
+		buf2 = buf + (i - line1) * ncols
+		call amovr (Memr[buf1], Memr[buf2], ncols)
+	    }
+	} else if (line2 > llast2) {
+	    do i = line1, line2 {
+		if (i < llast2)
+		    buf1 = buf + (i - llast1) * ncols
+		else
+		    buf1 = imgs2r (im, col1, col2, i, i)
+		buf2 = buf + (i - line1) * ncols
+		call amovr (Memr[buf1], Memr[buf2], ncols)
+	    }
+	}
+
+	llast1 = line1
+	llast2 = line2
+	nclast = ncols
+	nllast = nlines
+end
diff --git a/pkg/images/immatch/src/xregister/rgxplot.x b/pkg/images/immatch/src/xregister/rgxplot.x
new file mode 100644
index 00000000..8b347ab5
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxplot.x
@@ -0,0 +1,317 @@
+include <imhdr.h>
+include <gset.h>
+
+# RG_XPLINE -- Plot a line of reference and input image.
+
+procedure rg_xpline (gd, imr, im, nliner, xshift, yshift)
+
+pointer	gd		#I pointer to the graphics stream
+pointer	imr		#I pointer to the reference image
+pointer	im		#I pointer to the image
+int	nliner		#I the reference line
+int	xshift		#I x shift
+int	yshift		#I y shift
+
+int	i, rncols, rnlines, incols, inlines
+pointer	sp, title, xr, xi, ptrr, ptri
+real	ymin, ymax, tymin, tymax
+int	strlen()
+pointer	imgl1r(), imgl2r()
+
+begin
+	# Return if no graphics stream.
+	if (gd == NULL)
+	    return
+
+	# Check for valid line number.
+	rncols = IM_LEN(imr,1)
+	if (IM_NDIM(imr) == 1)
+	    rnlines = 1
+	else
+	    rnlines = IM_LEN(imr,2)
+	incols = IM_LEN(im,1)
+	if (IM_NDIM(im) == 1)
+	    inlines = 1
+	else
+	    inlines = IM_LEN(im,2)
+	if ((nliner < 1) || (nliner > rnlines))
+	    return
+	if (((nliner + yshift) < 1) || ((nliner + yshift) > inlines)) 
+	    return
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (title, SZ_LINE, TY_CHAR)
+	call salloc (xr, rncols, TY_REAL)
+	call salloc (xi, rncols, TY_REAL)
+
+	# Initialize the x data data.
+	do i = 1, rncols {
+	    Memr[xr+i-1] = i
+	    Memr[xi+i-1] = i - xshift
+	}
+
+	# Initalize the y data.
+	if (IM_NDIM(imr) == 1)
+	    ptrr = imgl1r (imr)
+	else
+	    ptrr = imgl2r (imr, nliner)
+	if (IM_NDIM(im) == 1)
+	    ptri = imgl1r (im)
+	else
+	    ptri = imgl2r (im, nliner + yshift)
+	call alimr (Memr[ptrr], rncols, ymin, ymax)
+	call alimr (Memr[ptri], incols, tymin, tymax)
+	ymin = min (ymin, tymin)
+	ymax = max (ymax, tymax)
+
+	# Construct the title.
+	call sprintf (Memc[title], SZ_LINE,
+	    "Refimage: %s  Image: %s\n")
+	    call pargstr (IM_HDRFILE(imr))
+	    call pargstr (IM_HDRFILE(im))
+	call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+	    "Refline (solid): %d  Inline (dashed): %d  Xlag: %d  Ylag: %d")
+	    call pargi (nliner)
+	    call pargi (nliner + yshift)
+	    call pargi (xshift)
+	    call pargi (yshift)
+
+	# Set up the axes labels and window.
+	call gclear (gd)
+	call gswind (gd, 1.0, real(rncols), ymin, ymax)
+	call glabax (gd, Memc[title], "Column Number", "Counts")
+
+	# Plot the two lines.
+	call gseti (gd, G_PLTYPE, GL_SOLID)
+	call gpline (gd, Memr[xr], Memr[ptrr], rncols)
+	call gseti (gd, G_PLTYPE, GL_DASHED)
+	call gpline (gd, Memr[xi], Memr[ptri], incols)
+	call gflush (gd)
+
+	call sfree (sp)
+end
+
+
+# RG_XPCOL -- Plot a column in the reference and input image.
+
+procedure rg_xpcol (gd, imr, im, ncolr, xshift, yshift)
+
+pointer	gd		#I pointer to the graphics stream
+pointer	imr		#I pointer to the reference image
+pointer	im		#I pointer to the image
+int	ncolr		#I the line number
+int	xshift		#I xshift to be applied
+int	yshift		#I yshift to be applied
+
+int	i, rncols, rnlines, incols, inlines
+pointer	sp, title, xr, xi, ptrr, ptri
+real	ymin, ymax, tymin, tymax
+int	strlen()
+pointer	imgs1r(), imgs2r()
+
+begin
+	# Return if no graphics stream.
+	if (gd == NULL)
+	    return
+
+	# Check for valid column number.
+	rncols = IM_LEN(imr,1)
+	if (IM_NDIM(imr) == 1)
+	    rnlines = 1
+	else
+	    rnlines = IM_LEN(imr,2)
+	incols = IM_LEN(im,1)
+	if (IM_NDIM(im) == 1)
+	    inlines = 1
+	else
+	    inlines = IM_LEN(im,2)
+	if ((ncolr < 1) || (ncolr > rncols))
+	    return
+	if (((ncolr - xshift) < 1) || ((ncolr - xshift) > incols))
+	    return
+
+	# Allocate valid working space.
+	call smark (sp)
+	call salloc (title, SZ_LINE, TY_CHAR)
+	call salloc (xr, rnlines, TY_REAL)
+	call salloc (xi, inlines, TY_REAL)
+
+	# Initialize the data.
+	do i = 1, rnlines {
+	    Memr[xr+i-1] = i
+	    Memr[xi+i-1] = i - yshift
+	}
+	if (IM_NDIM(imr) == 1)
+	    ptrr = imgs1r (imr, ncolr, ncolr)
+	else
+	    ptrr = imgs2r (imr, ncolr, ncolr, 1, rnlines)
+	if (IM_NDIM(im) == 1)
+	    ptri = imgs1r (im, ncolr + xshift, ncolr + xshift)
+	else
+	    ptri = imgs2r (im, ncolr + xshift, ncolr + xshift, 1, inlines)
+	call alimr (Memr[ptrr], rnlines, ymin, ymax)
+	call alimr (Memr[ptri], inlines, tymin, tymax)
+	ymin = min (ymin, tymin)
+	ymax = max (ymax, tymax)
+
+	# Construct the title.
+	call sprintf (Memc[title], SZ_LINE, "Refimage: %s Image: %s\n")
+	    call pargstr (IM_HDRFILE(imr))
+	    call pargstr (IM_HDRFILE(im))
+	call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+	    "Refcol (solid): %d  Imcol (dashed): %d  Xlag: %d Ylag: %d")
+	    call pargi (ncolr)
+	    call pargi (ncolr + xshift)
+	    call pargi (xshift)
+	    call pargi (yshift)
+
+	# Set up the labels and the axes.
+	call gclear (gd)
+	call gswind (gd, 1.0, real (rnlines), ymin, ymax)
+	call glabax (gd, Memc[title], "Line Number", "Counts")
+
+	# Plot the profile.
+	call gseti (gd, G_PLTYPE, GL_SOLID)
+	call gpline (gd, Memr[xr], Memr[ptrr], rnlines)
+	call gseti (gd, G_PLTYPE, GL_DASHED)
+	call gpline (gd, Memr[xi], Memr[ptri], rnlines)
+	call gflush (gd)
+
+	call sfree (sp)
+end
+
+
+# RG_XCPLINE -- Plot a line of the 2D correlation function.
+
+procedure rg_xcpline (gd, title, data, nx, ny, nline)
+
+pointer	gd		#I pointer to the graphics stream
+char	title[ARB]	#I title for the plot
+real	data[nx,ARB]	#I the input data array	
+int	nx, ny		#I dimensions of the input data array
+int	nline		#I the line number
+
+int	i
+pointer	sp, str, x
+real	ymin, ymax
+
+begin
+	# Return if no graphics stream.
+	if (gd == NULL)
+	    return
+
+	# Check for valid line number.
+	if (nline < 1 || nline > ny)
+	    return
+
+	# Allocate some working space.
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (x, nx, TY_REAL)
+
+	# Initialize the data.
+	do i = 1, nx
+	    Memr[x+i-1] = i
+	call alimr (data[1,nline], nx, ymin, ymax)
+
+	# Set up the labels and the axes.
+	call gclear (gd)
+	call gswind (gd, 1.0, real (nx), ymin, ymax)
+	call glabax (gd, title, "X Lag", "X-Correlation Function")
+
+	# Plot the line profile.
+	call gseti (gd, G_PLTYPE, GL_SOLID)
+	call gpline (gd, Memr[x], data[1,nline], nx)
+	call gflush (gd)
+
+	call sfree (sp)
+end
+
+
+# RG_XCPCOL -- Plot a column of the cross-correlation function.
+
+procedure rg_xcpcol (gd, title, data, nx, ny, ncol)
+
+pointer	gd		#I pointer to the graphics stream
+char	title[ARB]	#I title of the column plot
+real	data[nx,ARB]	#I the input data array
+int	nx, ny		#I the dimensions of the input data array
+int	ncol		#I line number
+
+int	i
+pointer	sp, x, y
+real	ymin, ymax
+
+begin
+	# Return if no graphics stream.
+	if (gd == NULL)
+	    return
+
+	# Check for valid column number.
+	if (ncol < 1 || ncol > nx)
+	    return
+
+	# Initialize.
+	call smark (sp)
+	call salloc (x, ny, TY_REAL)
+	call salloc (y, ny, TY_REAL)
+
+	# Get the data to be plotted.
+	do i = 1, ny {
+	    Memr[x+i-1] = i
+	    Memr[y+i-1] = data[ncol,i]
+	}
+	call alimr (Memr[y], ny, ymin, ymax)
+
+	# Set up the labels and the axes.
+	call gclear (gd)
+	call gswind (gd, 1.0, real (ny), ymin, ymax)
+	call glabax (gd, title, "Y Lag", "X-Correlation Function") 
+
+	# Plot the profile.
+	call gseti (gd, G_PLTYPE, GL_SOLID)
+	call gpline (gd, Memr[x], Memr[y], ny)
+	call gflush (gd)
+
+	call sfree (sp)
+end
+
+
+# RG_XMKPEAK -- Procedure to mark the peak from a correlation function
+# contour plot.
+
+procedure rg_xmkpeak (gd, xwindow, ywindow, xshift, yshift)
+
+pointer	gd		#I pointer to the graphics stream
+int	xwindow		#I x dimension of correlation function
+int	ywindow		#I y dimension of correlation function
+real	xshift		#O x shift
+real	yshift		#O y shift
+
+int	wcs, key
+pointer	sp, cmd
+real	wx, wy
+int	clgcur()
+
+begin
+	if (gd == NULL)
+	    return
+
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	call printf ("Mark peak of the cross correlation function\n")
+	if (clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) == EOF)
+	    ;
+	if (wx < 1.0 || wx > real (xwindow) || wy < 1.0 || wy >
+	    real (ywindow)) {
+	    xshift = 0.0
+	    yshift = 0.0
+	} else {
+	    xshift = wx - (1 + xwindow) / 2
+	    yshift = wy - (1 + ywindow) / 2
+	}
+
+	call sfree (sp)
+end
diff --git a/pkg/images/immatch/src/xregister/rgxppars.x b/pkg/images/immatch/src/xregister/rgxppars.x
new file mode 100644
index 00000000..2dc6aafd
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxppars.x
@@ -0,0 +1,49 @@
+include "xregister.h"
+
+# RG_PXPARS -- Update the cross-correlation algorithm parameters.
+
+procedure rg_pxpars (xc)
+
+pointer	xc		#I pointer to the cross-correlation structure
+
+pointer	sp, str
+int	rg_xstati()
+real	rg_xstatr()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Define the regions.
+	call rg_xstats (xc, REGIONS, Memc[str], SZ_LINE)
+	call clpstr ("regions", Memc[str])
+	call clputi ("xlag", rg_xstati (xc, XLAG))
+	call clputi ("ylag", rg_xstati (xc, YLAG))
+	call clputi ("dxlag", rg_xstati (xc, DXLAG))
+	call clputi ("dylag", rg_xstati (xc, DYLAG))
+
+	# Store the background fitting parameters.
+	call rg_xstats (xc, BSTRING, Memc[str], SZ_LINE)
+	call clpstr ("background", Memc[str])
+	call clputi ("border", rg_xstati (xc, BORDER))
+	call clputr ("loreject", rg_xstatr (xc, LOREJECT))
+	call clputr ("hireject", rg_xstatr (xc, HIREJECT))
+	call clputr ("apodize", rg_xstatr (xc, APODIZE))
+	call rg_xstats (xc, FSTRING, Memc[str], SZ_LINE)
+	call clpstr ("filter", Memc[str])
+
+	# Store the cross-correlation parameters.
+	call rg_xstats (xc, CSTRING, Memc[str], SZ_LINE)
+	call clpstr ("correlation", Memc[str])
+	call clputi ("xwindow", rg_xstati (xc, XWINDOW))
+	call clputi ("ywindow", rg_xstati (xc, YWINDOW))
+
+	# Store the peak centering parameters.
+	call rg_xstats (xc, PSTRING, Memc[str], SZ_LINE)
+	call clpstr ("function", Memc[str])
+	call clputi ("xcbox", rg_xstati (xc, XCBOX))
+	call clputi ("ycbox", rg_xstati (xc, YCBOX))
+
+	call sfree (sp)
+end
diff --git a/pkg/images/immatch/src/xregister/rgxregions.x b/pkg/images/immatch/src/xregister/rgxregions.x
new file mode 100644
index 00000000..ed682f61
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxregions.x
@@ -0,0 +1,459 @@
+include <fset.h>
+include <ctype.h>
+include <imhdr.h>
+include "xregister.h"
+
+# RG_XREGIONS -- Decode the image sections into regions. If the sections string
+# is NULL then the regions list is initially empty and depending on the mode
+# of the task, XREGISTER will or will not complain.Otherwise the image
+# sections specified in the sections string or file are decoded into a
+# regions list.
+
+int procedure rg_xregions (list, im, xc, rp)
+
+int	list			#I pointer to the regions list
+pointer	im			#I pointer to the reference image
+pointer	xc			#I pointer to the cross-correlation structure
+int	rp			#I index of the current region
+
+int	fd, nregions
+pointer	sp, fname, regions
+int	rg_xgrid(), rg_xgregions(), rg_xrregions(), rg_xstati(), fntgfnb()
+int	open()
+errchk	fntgfnb(), open(), close()
+
+begin
+	call smark (sp)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+	call salloc (regions, SZ_LINE, TY_CHAR)
+
+	call rg_xstats (xc, REGIONS, Memc[regions], SZ_LINE)
+	if (rp < 1 || rp > MAX_NREGIONS || Memc[regions] == EOS) {
+	    nregions = 0
+	} else if (rg_xgrid (im, xc, rp, MAX_NREGIONS) > 0) {
+	    nregions = rg_xstati (xc, NREGIONS)
+	} else if (rg_xgregions (im, xc, rp, MAX_NREGIONS) > 0) {
+	    nregions = rg_xstati (xc, NREGIONS)
+	} else if (list != NULL) {
+	    iferr {
+		if (fntgfnb (list, Memc[fname], SZ_FNAME) != EOF) {
+	            fd = open (Memc[fname], READ_ONLY, TEXT_FILE)
+	            nregions= rg_xrregions (fd, im, xc, rp, MAX_NREGIONS)
+	            call close (fd)
+		}
+	    } then
+		nregions = 0
+	} else
+	    nregions = 0
+
+	call sfree (sp)
+
+	return (nregions)
+end
+
+
+# RG_XMKREGIONS -- Create a list of regions by marking image sections
+# on the image display.
+
+int procedure rg_xmkregions (im, xc, rp, max_nregions, regions, maxch)
+
+pointer	im			#I pointer to the reference image
+pointer	xc			#I pointer to the cross-correlation structure
+int	rp			#I index of the current region
+int	max_nregions		#I the maximum number of regions
+char	regions[ARB]		#O the output regions string
+int	maxch			#I maximum size of the output regions string
+
+int	op, nregions, wcs, key
+pointer	sp, region, section, cmd
+real	xll, yll, xur, yur
+int	rg_xstati(), clgcur(), gstrcpy()
+pointer	rg_xstatp()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (region, SZ_LINE, TY_CHAR)
+	call salloc (section, SZ_LINE, TY_CHAR)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	# Allocate the arrays to hold the regions information,
+	call rg_xrealloc (xc, max_nregions)
+
+	# Initialize.
+	nregions = min (rp-1, rg_xstati (xc, NREGIONS))
+	op = 1
+
+	# Mark the sections on the display.
+	while (nregions < max_nregions) {
+
+	    call printf ("Mark lower left corner of region %d [q to quit].\n")
+		call pargi (nregions + 1)
+	    if (clgcur ("icommands", xll, yll, wcs, key, Memc[cmd],
+	        SZ_LINE) == EOF)
+		break
+	    if (key == 'q')
+		break
+
+	    call printf ("Mark upper right corner of region %d [q to quit].\n")
+		call pargi (nregions + 1)
+	    if (clgcur ("icommands", xur, yur, wcs, key, Memc[cmd],
+	        SZ_LINE) == EOF)
+		break
+	    if (key == 'q')
+		break
+
+	    if (xll < 1.0 || xur > IM_LEN(im,1) || yll < 1.0 || yur >
+		IM_LEN(im,2))
+		break
+
+	    Memi[rg_xstatp(xc,RC1)+nregions] = nint (xll)
+	    Memi[rg_xstatp(xc,RC2)+nregions] = nint (xur)
+	    Memi[rg_xstatp(xc,RL1)+nregions] = nint (yll)
+	    Memi[rg_xstatp(xc,RL2)+nregions] = nint (yur)
+	    Memr[rg_xstatp(xc,RZERO)+nregions] = INDEFR
+	    Memr[rg_xstatp(xc,RXSLOPE)+nregions] = INDEFR
+	    Memr[rg_xstatp(xc,RYSLOPE)+nregions] = INDEFR
+	    Memr[rg_xstatp(xc,XSHIFTS)+nregions] = INDEFR
+	    Memr[rg_xstatp(xc,YSHIFTS)+nregions] = INDEFR
+	    nregions = nregions + 1
+
+	    # Write the first 9 regions into the regions string.
+	    call sprintf (Memc[cmd], SZ_LINE, "[%d:%d,%d:%d] ")
+		call pargi (nint (xll))
+		call pargi (nint (xur))
+		call pargi (nint (yll))
+		call pargi (nint (yur))
+	    op = op + gstrcpy (Memc[cmd], regions[op], maxch - op + 1)
+	}
+	call printf ("\n")
+
+	# Reallocate the correct amount of space.
+	call rg_xseti (xc, NREGIONS, nregions)
+	if (nregions > 0)
+	    call rg_xrealloc (xc, nregions)
+	else 
+	    call rg_xrfree (xc)
+
+	call sfree (sp)
+
+	return (nregions)
+end
+
+
+# RG_XGRID - Decode the regions from a grid specification.
+
+int procedure rg_xgrid (im, xc, rp, max_nregions)
+
+pointer	im			#I pointer to the reference image
+pointer	xc			#I pointer to the cross-correlation structure
+int	rp			#I index of the current region
+int	max_nregions		#I the maximum number of regions
+
+int	i, istart, iend, j, jstart, jend, ncols, nlines, nxsample, nysample
+int	nxcols, nylines, nregions
+pointer	sp, region, section
+int	rg_xstati(), nscan(), strcmp()
+pointer	rg_xstatp()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (region, SZ_LINE, TY_CHAR)
+	call salloc (section, SZ_LINE, TY_CHAR)
+
+	# Allocate the arrays to hold the regions information,
+	call rg_xrealloc (xc, max_nregions)
+
+	# Initialize.
+	call rg_xstats (xc, REGIONS, Memc[region], SZ_LINE)
+	ncols = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	nregions = min (rp - 1, rg_xstati (xc, NREGIONS))
+
+	# Decode the grid specification.
+	call sscan (Memc[region])
+	    call gargwrd (Memc[section], SZ_LINE)
+	    call gargi (nxsample)
+	    call gargi (nysample) 
+	if ((nscan() != 3) || (strcmp (Memc[section], "grid") != 0)) {
+	    call sfree (sp)
+	    return (nregions)
+	}
+
+	# Decode the regions.
+	if ((nxsample * nysample) > max_nregions) {
+	    nxsample = nint (sqrt (real (max_nregions) * real (ncols) /
+	        real (nlines)))
+	    nysample = real (max_nregions) / real (nxsample)
+	}
+	nxcols = ncols / nxsample
+	nylines = nlines / nysample
+	jstart = 1 + (nlines - nysample * nylines) / 2
+	jend = jstart + (nysample - 1) * nylines
+	do j = jstart, jend, nylines {
+	    istart = 1 + (ncols - nxsample * nxcols) / 2
+	    iend = istart + (nxsample - 1) * nxcols
+	    do i = istart, iend, nxcols {
+		Memi[rg_xstatp(xc,RC1)+nregions] = i 
+		Memi[rg_xstatp(xc,RC2)+nregions] = i + nxcols - 1
+		Memi[rg_xstatp(xc,RL1)+nregions] = j
+		Memi[rg_xstatp(xc,RL2)+nregions] = j + nylines - 1
+		Memr[rg_xstatp(xc,RZERO)+nregions] = INDEFR
+		Memr[rg_xstatp(xc,RXSLOPE)+nregions] = INDEFR
+		Memr[rg_xstatp(xc,RYSLOPE)+nregions] = INDEFR
+		Memr[rg_xstatp(xc,XSHIFTS)+nregions] = INDEFR
+		Memr[rg_xstatp(xc,YSHIFTS)+nregions] = INDEFR
+		nregions = nregions + 1
+	    }
+	}
+
+	call rg_xseti (xc, NREGIONS, nregions)
+	if (nregions > 0)
+	    call rg_xrealloc (xc, nregions)
+	else
+	    call rg_xrfree (xc)
+	call sfree (sp)
+
+	return (nregions)
+end
+
+
+# RG_XRREGIONS -- Read and decode the regions from a file.
+
+int procedure rg_xrregions (fd, im, xc, rp, max_nregions)
+
+int	fd			#I regions file descriptor
+pointer	im			#I pointer to the reference image
+pointer	xc			#I pointer to the cross-correlation structure
+int	rp			#I index of the current region
+int	max_nregions		#I the maximum number of regions
+
+int	ncols, nlines, nregions, x1, y1, x2, y2, step
+pointer	sp, line, section
+int	rg_xstati(), getline(), rg_xgsections()
+pointer	rg_xstatp()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (line, SZ_LINE, TY_CHAR)
+	call salloc (section, SZ_LINE, TY_CHAR)
+
+	# Allocate the arrays to hold the regions information,
+	call rg_xrealloc (xc, max_nregions)
+
+	# Initialize.
+	ncols = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	nregions = min (rp - 1, rg_xstati (xc, NREGIONS))
+
+	# Decode the regions string.
+	while ((getline (fd, Memc[line]) != EOF) &&  nregions < max_nregions) {
+	    call sscan (Memc[line])
+		call gargwrd (Memc[section], SZ_LINE)
+	    while ((Memc[section] != EOS) && (nregions < max_nregions)) {
+		if (rg_xgsections (Memc[section], x1, x2, step, y1, y2, step,
+		    ncols, nlines) == OK) {
+		    Memi[rg_xstatp(xc,RC1)+nregions] = x1
+		    Memi[rg_xstatp(xc,RC2)+nregions] = x2
+		    Memi[rg_xstatp(xc,RL1)+nregions] = y1
+		    Memi[rg_xstatp(xc,RL2)+nregions] = y2
+		    Memr[rg_xstatp(xc,RZERO)+nregions] = INDEFR
+		    Memr[rg_xstatp(xc,RXSLOPE)+nregions] = INDEFR
+		    Memr[rg_xstatp(xc,RYSLOPE)+nregions] = INDEFR
+		    Memr[rg_xstatp(xc,XSHIFTS)+nregions] = INDEFR
+		    Memr[rg_xstatp(xc,YSHIFTS)+nregions] = INDEFR
+		    nregions = nregions + 1
+		}
+		call gargwrd (Memc[section], SZ_LINE)
+	    }
+	}
+
+	# Reallocate the correct amount of space.
+	call rg_xseti (xc, NREGIONS, nregions)
+	if (nregions > 0)
+	    call rg_xrealloc (xc, nregions)
+	else
+	    call rg_xrfree (xc)
+
+	call sfree (sp)
+
+	return (nregions)
+end
+
+
+# RG_XGREGIONS -- Decode a list of regions from a string containing
+# a list of sections.
+
+int procedure rg_xgregions (im, xc, rp, max_nregions)
+
+pointer	im			#I pointer to the reference image
+pointer	xc			#I pointer to cross-correlation structure
+int	rp			#I the index of the current region
+int	max_nregions		#I the maximum number of regions
+
+int	ncols, nlines, nregions, x1, x2, y1, y2, step
+pointer	sp, section, region
+int	rg_xstati(), rg_xgsections()
+pointer	rg_xstatp()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (region, SZ_LINE, TY_CHAR)
+	call salloc (section, SZ_LINE, TY_CHAR)
+
+	# Allocate the arrays to hold the regions information.
+	call rg_xrealloc (xc, max_nregions)
+
+	# Initialize.
+	call rg_xstats (xc, REGIONS, Memc[region], SZ_LINE)
+	ncols = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	nregions = min (rp - 1, rg_xstati (xc, NREGIONS))
+
+	# Decode the sections
+	call sscan (Memc[region])
+	    call gargwrd (Memc[section], SZ_LINE)
+	while ((Memc[section] != EOS) && (nregions < max_nregions)) {
+	    if (rg_xgsections (Memc[section], x1, x2, step, y1, y2, step,
+		ncols, nlines) == OK) {
+		Memi[rg_xstatp(xc,RC1)+nregions] = x1
+		Memi[rg_xstatp(xc,RC2)+nregions] = x2
+		Memi[rg_xstatp(xc,RL1)+nregions] = y1
+		Memi[rg_xstatp(xc,RL2)+nregions] = y2
+		Memr[rg_xstatp(xc,RZERO)+nregions] = INDEFR
+		Memr[rg_xstatp(xc,RXSLOPE)+nregions] = INDEFR
+		Memr[rg_xstatp(xc,RYSLOPE)+nregions] = INDEFR
+		Memr[rg_xstatp(xc,XSHIFTS)+nregions] = INDEFR
+		Memr[rg_xstatp(xc,YSHIFTS)+nregions] = INDEFR
+		nregions = nregions + 1
+	    }
+	    call gargwrd (Memc[section], SZ_LINE)
+	}
+
+
+	# Reallocate the correct amount of space.
+	call rg_xseti (xc, NREGIONS, nregions)
+	if (nregions > 0)
+	    call rg_xrealloc (xc, nregions)
+	else
+	    call rg_xrfree (xc)
+
+	call sfree (sp)
+
+	return (nregions)
+end
+
+
+# RG_XGSECTIONS -- Decode an image section into column and line limits
+# and a step size. Sections which describe the whole image are decoded into
+# a block ncols * nlines long.
+
+int procedure rg_xgsections (section, x1, x2, xstep, y1, y2, ystep, ncols,
+	nlines)
+
+char	section[ARB]		#I the input section string
+int	x1, x2			#O the output column section limits
+int	xstep			#O the output column step size
+int	y1, y2			#O the output line section limits
+int	ystep			#O the output line step size
+int	ncols, nlines		#I the maximum number of lines and columns
+
+int	ip
+int	rg_xgdim()
+
+begin
+	ip = 1
+	if (rg_xgdim (section, ip, x1, x2, xstep, ncols) == ERR)
+	    return (ERR)
+	if (rg_xgdim (section, ip, y1, y2, ystep, nlines) == ERR)
+	    return (ERR)
+
+	return (OK)
+end
+
+
+# RG_XGDIM -- Decode a single subscript expression to produce the
+# range of values for that subscript (X1:X2), and the sampling step size, STEP.
+# Note that X1 may be less than, greater than, or equal to X2, and STEP may
+# be a positive or negative nonzero integer.  Various shorthand notations are
+# permitted, as is embedded whitespace.
+
+int procedure rg_xgdim (section, ip, x1, x2, step, limit)
+
+char	section[ARB]		#I the input image section
+int	ip			#I/O pointer to the position in section string
+int	x1			#O first limit of dimension
+int	x2			#O second limit of dimension
+int	step			#O step size of dimension
+int	limit			#I maximum size of dimension
+
+int	temp
+int	ctoi()
+
+begin
+	x1 = 1
+	x2 = limit
+	step = 1
+
+	while (IS_WHITE(section[ip]))
+	    ip = ip + 1
+
+	if (section[ip] =='[')
+	    ip = ip + 1
+
+	while (IS_WHITE(section[ip]))
+	    ip = ip + 1
+
+	# Get X1, X2.
+	if (ctoi (section, ip, temp) > 0) {			# [x1
+	    x1 = max (1, min (temp, limit))
+	    if (section[ip] == ':') {	
+		ip = ip + 1
+		if (ctoi (section, ip, temp) == 0)		# [x1:x2
+		    return (ERR)
+		x2 = max (1, min (temp, limit))
+	    } else
+		x2 = x1
+
+	} else if (section[ip] == '-') {
+	    x1 = limit
+	    x2 = 1
+	    ip = ip + 1
+	    if (section[ip] == '*')
+		ip = ip + 1
+
+	} else if (section[ip] == '*')				# [*
+	    ip = ip + 1
+
+	while (IS_WHITE(section[ip]))
+	    ip = ip + 1
+
+	# Get sample step size, if give.
+	if (section[ip] == ':') {				# ..:step
+	    ip = ip + 1
+	    if (ctoi (section, ip, step) == 0)
+		return (ERR)
+	    else if (step == 0)
+		return (ERR)
+	}
+
+	# Allow notation such as "-*:5", (or even "-:5") where the step
+	# is obviously supposed to be negative.
+
+	if (x1 > x2 && step > 0)
+	    step = -step
+	
+	while (IS_WHITE(section[ip]))
+	    ip = ip + 1
+
+	if (section[ip] == ',') {
+	    ip = ip + 1
+	    return (OK)
+	} else if (section[ip] == ']')
+	    return (OK)
+	else
+	    return (ERR)
+end
diff --git a/pkg/images/immatch/src/xregister/rgxshow.x b/pkg/images/immatch/src/xregister/rgxshow.x
new file mode 100644
index 00000000..3a746d9c
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxshow.x
@@ -0,0 +1,172 @@
+include "xregister.h"
+
+# RG_XSHOW -- Show the XREGISTER parameters.
+
+procedure rg_xshow (xc)
+
+pointer	xc		#I pointer to the main xregister structure
+
+begin
+	call rg_xnshow (xc)
+	call printf ("\n")
+	call rg_xbshow (xc)
+	call printf ("\n")
+	call rg_xxshow (xc)
+	call printf ("\n")
+	call rg_xpshow (xc)
+end
+
+
+# RG_XNSHOW -- Show the input/output data XREGISTER parameters.
+
+procedure rg_xnshow (xc)
+
+pointer	xc	#I pointer to the main xregister structure
+
+pointer	sp, str
+int	rg_xstati()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Set the object characteristics.
+	call printf ("\nInput/output data\n")
+	call rg_xstats (xc, IMAGE, Memc[str], SZ_FNAME)
+	call printf ("    %s: %s\n")
+	    call pargstr (KY_IMAGE)
+	    call pargstr (Memc[str])
+	call rg_xstats (xc, REFIMAGE, Memc[str], SZ_FNAME)
+	call printf ("    %s: %s\n")
+	    call pargstr (KY_REFIMAGE)
+	    call pargstr (Memc[str])
+	call rg_xstats (xc, REGIONS, Memc[str], SZ_FNAME)
+	call printf ("    %s: %s\n")
+	    call pargstr (KY_REGIONS)
+	    call pargstr (Memc[str])
+	call printf ("    %s = %d    %s = %d\n")
+	    call pargstr (KY_XLAG)
+	    call pargi (rg_xstati (xc, XLAG))
+	    call pargstr (KY_YLAG)
+	    call pargi (rg_xstati (xc, YLAG))
+	call printf ("    %s = %d    %s = %d\n")
+	    call pargstr (KY_DXLAG)
+	    call pargi (rg_xstati (xc, DXLAG))
+	    call pargstr (KY_DYLAG)
+	    call pargi (rg_xstati (xc, DYLAG))
+	call rg_xstats (xc, DATABASE, Memc[str], SZ_FNAME)
+	call printf ("    %s: %s\n")
+	    call pargstr (KY_DATABASE)
+	    call pargstr (Memc[str])
+	call rg_xstats (xc, RECORD, Memc[str], SZ_FNAME)
+	call printf ("    %s: %s\n")
+	    call pargstr (KY_RECORD)
+	    call pargstr (Memc[str])
+	call rg_xstats (xc, REFFILE, Memc[str], SZ_FNAME)
+	call printf ("    %s: %s\n")
+	    call pargstr (KY_REFFILE)
+	    call pargstr (Memc[str])
+	call rg_xstats (xc, OUTIMAGE, Memc[str], SZ_FNAME)
+	call printf ("    %s: %s\n")
+	    call pargstr (KY_OUTIMAGE)
+	    call pargstr (Memc[str])
+
+	call sfree (sp)
+end
+
+
+# RG_XBSHOW -- Show the background fitting parameters.
+
+procedure rg_xbshow (xc)
+
+pointer	xc		#I pointer to the main xregister structure
+
+int	back
+pointer	sp, str
+int	rg_xstati()
+real	rg_xstatr()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	back = rg_xstati (xc, BACKGRD)
+	call printf ("Background fitting parameters:\n")
+	call rg_xstats (xc, BSTRING, Memc[str], SZ_LINE)
+	call printf ("    %s:  %s\n")
+	    call pargstr (KY_BACKGROUND)
+	    call pargstr (Memc[str])
+	call printf ("    %s = %d\n")
+	    call pargstr (KY_BORDER)
+	    call pargi (rg_xstati (xc, BORDER))
+	call printf ("    %s = %g   %s = %g\n")
+	    call pargstr (KY_LOREJECT)
+	    call pargr (rg_xstatr (xc, LOREJECT))
+	    call pargstr (KY_HIREJECT)
+	    call pargr (rg_xstatr (xc, HIREJECT))
+	call printf ("    %s = %g\n")
+	    call pargstr (KY_APODIZE)
+	    call pargr (rg_xstatr (xc, APODIZE))
+	call rg_xstats (xc, FSTRING, Memc[str], SZ_LINE)
+	call printf ("    %s:  %s\n")
+	    call pargstr (KY_FILTER)
+	    call pargstr (Memc[str])
+
+	call sfree (sp)
+end
+
+
+# RG_XXSHOW -- Show the cross-correlation function parameters.
+
+procedure rg_xxshow (xc)
+
+pointer	xc		#I pointer to the main xregister structure
+
+pointer	sp, str
+int	rg_xstati()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	call printf ("Cross correlation function:\n")
+	call rg_xstats (xc, CSTRING, Memc[str], SZ_LINE)
+	call printf ("    %s:  %s\n")
+	    call pargstr (KY_CORRELATION)
+	    call pargstr (Memc[str])
+	call printf ("    %s = %d    %s = %d\n")
+	    call pargstr (KY_XWINDOW)
+	    call pargi (rg_xstati (xc, XWINDOW))
+	    call pargstr (KY_YWINDOW)
+	    call pargi (rg_xstati (xc, YWINDOW))
+
+	call sfree (sp)
+end
+
+
+# RG_XPSHOW -- Show the peak centering parameters.
+
+procedure rg_xpshow (xc)
+
+pointer	xc		#I pointer to the main xregister structure
+
+pointer	sp, str
+int	rg_xstati()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	call printf ("Peak centering parameters:\n")
+	call rg_xstats (xc, PSTRING, Memc[str], SZ_LINE)
+	call printf ("    %s:  %s\n")
+	    call pargstr (KY_PEAKCENTER)
+	    call pargstr (Memc[str])
+	call printf ("    %s = %d    %s = %d\n")
+	    call pargstr (KY_XCBOX)
+	    call pargi (rg_xstati (xc, XCBOX))
+	    call pargstr (KY_YCBOX)
+	    call pargi (rg_xstati (xc, YCBOX))
+
+	call sfree (sp)
+end
diff --git a/pkg/images/immatch/src/xregister/rgxtools.x b/pkg/images/immatch/src/xregister/rgxtools.x
new file mode 100644
index 00000000..e1fb921e
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxtools.x
@@ -0,0 +1,685 @@
+include "xregister.h"
+
+#  RG_XINIT -- Initialize the cross-correlation code fitting structure. 
+
+procedure rg_xinit (xc, cfunc)
+
+pointer	xc		#O pointer to the cross-correlation structure
+int	cfunc		#I the input cross-correlation function
+
+begin
+	call malloc (xc, LEN_XCSTRUCT, TY_STRUCT)
+
+	# Initialize the regions pointers.
+	XC_RC1(xc) = NULL
+	XC_RC2(xc) = NULL
+	XC_RL1(xc) = NULL
+	XC_RL2(xc) = NULL
+	XC_RZERO(xc) = NULL
+	XC_RXSLOPE(xc) = NULL
+	XC_RYSLOPE(xc) = NULL
+	XC_XSHIFTS(xc) = NULL
+	XC_YSHIFTS(xc) = NULL
+	XC_TXSHIFT(xc) = 0.0
+	XC_TYSHIFT(xc) = 0.0
+	XC_NREGIONS(xc) = 0
+	XC_CREGION(xc) = 1
+
+	# Set up transformation parameters.
+	XC_NREFPTS(xc) = 0
+	call malloc (XC_XREF(xc), MAX_NREF, TY_REAL)
+	call malloc (XC_YREF(xc), MAX_NREF, TY_REAL)
+	call malloc (XC_TRANSFORM(xc), MAX_NTRANSFORM, TY_REAL)
+
+	# Initialize the region offsets
+	XC_IXLAG(xc) = DEF_IXLAG
+	XC_IYLAG(xc) = DEF_IYLAG
+	XC_XLAG(xc) = DEF_IXLAG
+	XC_YLAG(xc) = DEF_IYLAG
+	XC_DXLAG(xc) = DEF_DXLAG
+	XC_DYLAG(xc) = DEF_DYLAG
+		
+	# Define the background fitting parameters.
+	XC_BACKGRD(xc) = XC_BNONE
+	call strcpy ("none", XC_BSTRING(xc), SZ_FNAME)
+	XC_BVALUER(xc) = 0.0
+	XC_BVALUE(xc) = 0.0
+	XC_BORDER(xc) = DEF_BORDER
+	XC_LOREJECT(xc) = DEF_LOREJECT
+	XC_HIREJECT(xc) = DEF_HIREJECT
+	XC_APODIZE(xc) = 0.0
+	XC_FILTER(xc) = XC_FNONE
+	call strcpy ("none", XC_FSTRING(xc), SZ_FNAME)
+
+	# Get the correlation parameters.
+	XC_CFUNC(xc) = cfunc
+	switch (cfunc) {
+	case XC_DISCRETE:
+	    call strcpy ("discrete", XC_CSTRING(xc), SZ_FNAME)
+	case XC_FOURIER:
+	    call strcpy ("fourier", XC_CSTRING(xc), SZ_FNAME)
+	case XC_FILE:
+	    call strcpy ("file", XC_CSTRING(xc), SZ_FNAME)
+	case XC_DIFFERENCE:
+	    call strcpy ("difference", XC_CSTRING(xc), SZ_FNAME)
+	default:
+	    call strcpy ("unknown", XC_CSTRING(xc), SZ_FNAME)
+	}
+	XC_XWINDOW(xc) = DEF_XWINDOW
+	XC_YWINDOW(xc) = DEF_YWINDOW
+	XC_XCOR(xc) = NULL
+
+	# Define the peak fitting function.
+	XC_PFUNC(xc) = DEF_PFUNC
+	call sprintf (XC_PSTRING(xc), SZ_FNAME, "%s")
+	    call pargstr ("centroid")
+	XC_XCBOX(xc) = DEF_XCBOX
+	XC_YCBOX(xc) = DEF_YCBOX
+
+	# Initialize the strings.
+	XC_IMAGE(xc) = EOS
+	XC_REFIMAGE(xc) = EOS
+	XC_REGIONS(xc) = EOS
+	XC_DATABASE(xc) = EOS
+	XC_OUTIMAGE(xc) = EOS
+	XC_REFFILE(xc) = EOS
+	XC_RECORD(xc) = EOS
+
+	# Initialize the buffers.
+	call rg_xrinit (xc)
+
+end
+
+
+#  RG_XRINIT -- Initialize the regions definition portion of the 
+# cross correlation code fitting structure. 
+
+procedure rg_xrinit (xc)
+
+pointer	xc		#I pointer to crosscor structure
+
+begin
+	call rg_xrfree (xc)
+
+	XC_NREGIONS(xc) = 0
+	XC_CREGION(xc) = 1
+
+	call malloc (XC_RC1(xc), MAX_NREGIONS, TY_INT)
+	call malloc (XC_RC2(xc), MAX_NREGIONS, TY_INT)
+	call malloc (XC_RL1(xc), MAX_NREGIONS, TY_INT)
+	call malloc (XC_RL2(xc), MAX_NREGIONS, TY_INT)
+	call malloc (XC_RZERO(xc), MAX_NREGIONS, TY_REAL)
+	call malloc (XC_RXSLOPE(xc), MAX_NREGIONS, TY_REAL)
+	call malloc (XC_RYSLOPE(xc), MAX_NREGIONS, TY_REAL)
+	call malloc (XC_XSHIFTS(xc), MAX_NREGIONS, TY_REAL)
+	call malloc (XC_YSHIFTS(xc), MAX_NREGIONS, TY_REAL)
+
+	call amovki (INDEFI, Memi[XC_RC1(xc)], MAX_NREGIONS)
+	call amovki (INDEFI, Memi[XC_RC2(xc)], MAX_NREGIONS)
+	call amovki (INDEFI, Memi[XC_RL1(xc)], MAX_NREGIONS)
+	call amovki (INDEFI, Memi[XC_RL2(xc)], MAX_NREGIONS)
+	call amovkr (INDEFR, Memr[XC_RZERO(xc)], MAX_NREGIONS)
+	call amovkr (INDEFR, Memr[XC_RXSLOPE(xc)], MAX_NREGIONS)
+	call amovkr (INDEFR, Memr[XC_RYSLOPE(xc)], MAX_NREGIONS)
+	call amovkr (INDEFR, Memr[XC_XSHIFTS(xc)], MAX_NREGIONS)
+	call amovkr (INDEFR, Memr[XC_YSHIFTS(xc)], MAX_NREGIONS)
+
+	XC_TXSHIFT(xc) = 0.0
+	XC_TYSHIFT(xc) = 0.0
+end
+
+
+# RG_XCINDEFR -- Re-initialize the background and answers regions portion of
+# the cross-correlation fitting structure
+
+procedure rg_xcindefr (xc, creg)
+
+pointer	xc		#I pointer to the cross-correlation structure
+int     creg            #I the current region
+
+int	nregions
+int	rg_xstati()
+
+begin
+	nregions = rg_xstati (xc, NREGIONS)
+	if (creg < 1 || creg > nregions)
+	    return
+
+	if (nregions > 0) {
+	    Memr[XC_RZERO(xc)+creg-1] = INDEFR
+	    Memr[XC_RXSLOPE(xc)+creg-1] = INDEFR
+	    Memr[XC_RYSLOPE(xc)+creg-1] = INDEFR
+	    Memr[XC_XSHIFTS(xc)+creg-1] = INDEFR
+	    Memr[XC_YSHIFTS(xc)+creg-1] = INDEFR
+	}
+
+	XC_TXSHIFT(xc) = 0.0
+	XC_TYSHIFT(xc) = 0.0
+end
+
+
+# RG_XINDEFR -- Re-initialize the background and answers regions portion of
+# the cross-correlation fitting structure for all regions and reset the
+# current region to 1.
+
+procedure rg_xindefr (xc)
+
+pointer xc              #I pointer to the cross-correlation structure
+
+int     nregions
+int     rg_xstati()
+
+begin
+        nregions = rg_xstati (xc, NREGIONS)
+
+        if (nregions > 0) {
+            call amovkr (INDEFR, Memr[XC_RZERO(xc)], nregions)
+            call amovkr (INDEFR, Memr[XC_RXSLOPE(xc)], nregions)
+            call amovkr (INDEFR, Memr[XC_RYSLOPE(xc)], nregions)
+            call amovkr (INDEFR, Memr[XC_XSHIFTS(xc)], nregions)
+            call amovkr (INDEFR, Memr[XC_YSHIFTS(xc)], nregions)
+        }
+
+        XC_CREGION(xc) = 1
+        XC_TXSHIFT(xc) = 0.0
+        XC_TYSHIFT(xc) = 0.0
+end
+
+
+#  RG_XREALLOC -- Reallocate the regions bufffers and initialize if necessary.
+
+procedure rg_xrealloc (xc, nregions)
+
+pointer	xc		#I pointer to crosscor structure
+int	nregions	#I number of regions
+
+int	nr
+int	rg_xstati()
+
+begin
+	nr = rg_xstati (xc, NREGIONS)
+
+	call realloc (XC_RC1(xc), nregions, TY_INT)
+	call realloc (XC_RC2(xc), nregions, TY_INT)
+	call realloc (XC_RL1(xc), nregions, TY_INT)
+	call realloc (XC_RL2(xc), nregions, TY_INT)
+	call realloc (XC_RZERO(xc), nregions, TY_REAL)
+	call realloc (XC_RXSLOPE(xc), nregions, TY_REAL)
+	call realloc (XC_RYSLOPE(xc), nregions, TY_REAL)
+	call realloc (XC_XSHIFTS(xc), nregions, TY_REAL)
+	call realloc (XC_YSHIFTS(xc), nregions, TY_REAL)
+
+	call amovki (INDEFI, Memi[XC_RC1(xc)+nr], nregions - nr)
+	call amovki (INDEFI, Memi[XC_RC2(xc)+nr], nregions - nr)
+	call amovki (INDEFI, Memi[XC_RL1(xc)+nr], nregions - nr)
+	call amovki (INDEFI, Memi[XC_RL2(xc)+nr], nregions - nr)
+	call amovkr (INDEFR, Memr[XC_RZERO(xc)+nr], nregions - nr)
+	call amovkr (INDEFR, Memr[XC_RXSLOPE(xc)+nr], nregions - nr)
+	call amovkr (INDEFR, Memr[XC_RYSLOPE(xc)+nr], nregions - nr)
+	call amovkr (INDEFR, Memr[XC_XSHIFTS(xc)+nr], nregions - nr)
+	call amovkr (INDEFR, Memr[XC_YSHIFTS(xc)+nr], nregions - nr)
+end
+
+
+# RG_XFREE -- Free the cross-correlation fitting structure.
+
+procedure rg_xfree (xc)
+
+pointer	xc		#I pointer to the cross-correlation structure
+
+begin
+	# Free the region descriptors.
+	call rg_xrfree (xc)
+
+	# Free the transformation descriptors.
+	if (XC_XREF(xc) != NULL)
+	    call mfree (XC_XREF(xc), TY_REAL)
+	if (XC_YREF(xc) != NULL)
+	    call mfree (XC_YREF(xc), TY_REAL)
+	if (XC_TRANSFORM(xc) != NULL)
+	    call mfree (XC_TRANSFORM(xc), TY_REAL)
+
+	# Free the correlation function.
+	if (XC_XCOR(xc) != NULL)
+	    call mfree (XC_XCOR(xc), TY_REAL)
+
+	call mfree (xc, TY_STRUCT)
+end
+
+
+# RG_XRFREE -- Free the regions portion of the cross-correlation structure.
+
+procedure rg_xrfree (xc)
+
+pointer	xc		#I pointer to the cross-correlation structure
+
+begin
+	call rg_xseti (xc, NREGIONS, 0)
+	if (XC_RC1(xc) != NULL)
+	    call mfree (XC_RC1(xc), TY_INT)
+	XC_RC1(xc) = NULL
+	if (XC_RC2(xc) != NULL)
+	    call mfree (XC_RC2(xc), TY_INT)
+	XC_RC2(xc) = NULL
+	if (XC_RL1(xc) != NULL)
+	    call mfree (XC_RL1(xc), TY_INT)
+	XC_RL1(xc) = NULL
+	if (XC_RL2(xc) != NULL)
+	    call mfree (XC_RL2(xc), TY_INT)
+	XC_RL2(xc) = NULL
+	if (XC_RZERO(xc) != NULL)
+	    call mfree (XC_RZERO(xc), TY_REAL)
+	XC_RZERO(xc) = NULL
+	if (XC_RXSLOPE(xc) != NULL)
+	    call mfree (XC_RXSLOPE(xc), TY_REAL)
+	XC_RXSLOPE(xc) = NULL
+	if (XC_RYSLOPE(xc) != NULL)
+	    call mfree (XC_RYSLOPE(xc), TY_REAL)
+	XC_RYSLOPE(xc) = NULL
+	if (XC_XSHIFTS(xc) != NULL)
+	    call mfree (XC_XSHIFTS(xc), TY_REAL)
+	XC_XSHIFTS(xc) = NULL
+	if (XC_YSHIFTS(xc) != NULL)
+	    call mfree (XC_YSHIFTS(xc), TY_REAL)
+	XC_YSHIFTS(xc) = NULL
+end
+
+
+# RG_XSTATI -- Fetch the value of a cross-correlation fitting structure
+# integer parameter.
+
+int procedure rg_xstati (xc, param)
+
+pointer	xc		#I pointer to the cross-correlation fitting structure
+int	param		#I parameter to be fetched
+
+begin
+	switch (param) {
+	case CFUNC:
+	    return (XC_CFUNC(xc))
+	case IXLAG:
+	    return (XC_IXLAG(xc))
+	case IYLAG:
+	    return (XC_IYLAG(xc))
+	case XLAG:
+	    return (XC_XLAG(xc))
+	case YLAG:
+	    return (XC_YLAG(xc))
+	case DXLAG:
+	    return (XC_DXLAG(xc))
+	case DYLAG:
+	    return (XC_DYLAG(xc))
+	case XWINDOW:
+	    return (XC_XWINDOW(xc))
+	case YWINDOW:
+	    return (XC_YWINDOW(xc))
+	case CREGION:
+	    return (XC_CREGION(xc))
+	case NREGIONS:
+	    return (XC_NREGIONS(xc))
+	case BACKGRD:
+	    return (XC_BACKGRD(xc))
+	case BORDER:
+	    return (XC_BORDER(xc))
+	case FILTER:
+	    return (XC_FILTER(xc))
+	case XCBOX:
+	    return (XC_XCBOX(xc))
+	case YCBOX:
+	    return (XC_YCBOX(xc))
+	case PFUNC:
+	    return (XC_PFUNC(xc))
+	case NREFPTS:
+	    return (XC_NREFPTS(xc))
+	default:
+	    call error (0, "RG_XSTATI: Undefined integer parameter.")
+	}
+end
+
+
+# RG_XSTATP -- Fetch the value of a pointer parameter.
+
+pointer procedure rg_xstatp (xc, param)
+
+pointer	xc		#I pointer to the cross-correlation structure
+int	param		#I parameter to be fetched
+
+begin
+	switch (param) {
+	case RC1:
+	    return (XC_RC1(xc))
+	case RC2:
+	    return (XC_RC2(xc))
+	case RL1:
+	    return (XC_RL1(xc))
+	case RL2:
+	    return (XC_RL2(xc))
+	case RZERO:
+	    return (XC_RZERO(xc))
+	case RXSLOPE:
+	    return (XC_RXSLOPE(xc))
+	case RYSLOPE:
+	    return (XC_RYSLOPE(xc))
+	case XSHIFTS:
+	    return (XC_XSHIFTS(xc))
+	case YSHIFTS:
+	    return (XC_YSHIFTS(xc))
+	case XCOR:
+	    return (XC_XCOR(xc))
+	case XREF:
+	    return (XC_XREF(xc))
+	case YREF:
+	    return (XC_YREF(xc))
+#	case CORAPODIZE:
+#	    return (XC_CORAPODIZE(xc))
+	case TRANSFORM:
+	    return (XC_TRANSFORM(xc))
+	default:
+	    call error (0, "RG_XSTATP: Undefined pointer parameter.")
+	}
+end
+
+
+# RG_XSTATR -- Fetch the value of a real parameter.
+
+real procedure rg_xstatr (xc, param)
+
+pointer	xc		#I pointer to the cross-correlation structure
+int	param		#I parameter to be fetched
+
+begin
+	switch (param) {
+	case BVALUER:
+	    return (XC_BVALUER(xc))
+	case BVALUE:
+	    return (XC_BVALUE(xc))
+	case LOREJECT:
+	    return (XC_LOREJECT(xc))
+	case HIREJECT:
+	    return (XC_HIREJECT(xc))
+	case APODIZE:
+	    return (XC_APODIZE(xc))
+	case TXSHIFT:
+	    return (XC_TXSHIFT(xc))
+	case TYSHIFT:
+	    return (XC_TYSHIFT(xc))
+	default:
+	    call error (0, "RG_XSTATR: Undefined real parameter.")
+	}
+end
+
+
+# RG_XSTATS -- Fetch the value of a string parameter.
+
+procedure rg_xstats (xc, param, str, maxch)
+
+pointer	xc		#I pointer to the cross-correlation structure
+int	param		#I parameter to be fetched
+char	str[ARB]	#O output value of string parameter
+int	maxch		#I maximum number of characters in output string
+
+begin
+	switch (param) {
+	case BSTRING:
+	    call strcpy (XC_BSTRING(xc), str, maxch)
+	case FSTRING:
+	    call strcpy (XC_FSTRING(xc), str, maxch)
+	case CSTRING:
+	    call strcpy (XC_CSTRING(xc), str, maxch)
+	case PSTRING:
+	    call strcpy (XC_PSTRING(xc), str, maxch)
+	case REFIMAGE:
+	    call strcpy (XC_REFIMAGE(xc), str, maxch)
+	case IMAGE:
+	    call strcpy (XC_IMAGE(xc), str, maxch)
+	case OUTIMAGE:
+	    call strcpy (XC_OUTIMAGE(xc), str, maxch)
+	case REGIONS:
+	    call strcpy (XC_REGIONS(xc), str, maxch)
+	case DATABASE:
+	    call strcpy (XC_DATABASE(xc), str, maxch)
+	case RECORD:
+	    call strcpy (XC_RECORD(xc), str, maxch)
+	case REFFILE:
+	    call strcpy (XC_REFFILE(xc), str, maxch)
+	default:
+	    call error (0, "RG_XSTATS: Undefined string parameter.")
+	}
+end
+
+
+# RG_XSETI -- Set the value of an integer parameter.
+
+procedure rg_xseti (xc, param, value)
+
+pointer	xc		#I pointer to the cross-correlation structure
+int	param		#I parameter to be set
+int	value		#O value of the integer parameter
+
+begin
+	switch (param) {
+	case CFUNC:
+	    XC_CFUNC(xc) = value
+	    switch (value) {
+	    case XC_DISCRETE:
+		call strcpy ("discrete", XC_CSTRING(xc), SZ_FNAME)
+	    case XC_FOURIER:
+		call strcpy ("fourier", XC_CSTRING(xc), SZ_FNAME)
+	    case XC_FILE:
+		call strcpy ("file", XC_CSTRING(xc), SZ_FNAME)
+	    case XC_DIFFERENCE:
+		call strcpy ("difference", XC_CSTRING(xc), SZ_FNAME)
+	    default:
+		call strcpy ("unknown", XC_CSTRING(xc), SZ_FNAME)
+	    }
+	case IXLAG:
+	    XC_IXLAG(xc) = value
+	case IYLAG:
+	    XC_IYLAG(xc) = value
+	case XLAG:
+	    XC_XLAG(xc) = value
+	case YLAG:
+	    XC_YLAG(xc) = value
+	case DXLAG:
+	    XC_DXLAG(xc) = value
+	case DYLAG:
+	    XC_DYLAG(xc) = value
+	case XWINDOW:
+	    XC_XWINDOW(xc) = value
+	case YWINDOW:
+	    XC_YWINDOW(xc) = value
+	case BACKGRD:
+	    XC_BACKGRD(xc) = value
+	    switch (value) {
+	    case XC_BNONE:
+		call strcpy ("none", XC_BSTRING(xc), SZ_FNAME)
+	    case XC_MEAN:
+		call strcpy ("mean", XC_BSTRING(xc), SZ_FNAME)
+	    case XC_MEDIAN:
+		call strcpy ("median", XC_BSTRING(xc), SZ_FNAME)
+	    case XC_SLOPE:
+		call strcpy ("plane", XC_BSTRING(xc), SZ_FNAME)
+	    default:
+		call strcpy ("none", XC_BSTRING(xc), SZ_FNAME)
+	    }
+	case BORDER:
+	    XC_BORDER(xc) = value
+	case FILTER:
+	    XC_FILTER(xc) = value
+	    switch (value) {
+	    case XC_FNONE:
+		call strcpy ("none", XC_FSTRING(xc), SZ_FNAME)
+	    case XC_LAPLACE:
+		call strcpy ("laplace", XC_FSTRING(xc), SZ_FNAME)
+	    default:
+		call strcpy ("none", XC_FSTRING(xc), SZ_FNAME)
+	    }
+	case XCBOX:
+	    XC_XCBOX(xc) = value
+	case YCBOX:
+	    XC_YCBOX(xc) = value
+	case PFUNC:
+	    XC_PFUNC(xc) = value
+	    switch (value) {
+	    case XC_PNONE:
+		call strcpy ("none", XC_PSTRING(xc), SZ_FNAME)
+	    case XC_CENTROID:
+		call strcpy ("centroid", XC_PSTRING(xc), SZ_FNAME)
+	    case XC_PARABOLA:
+		call strcpy ("parabolic", XC_PSTRING(xc), SZ_FNAME)
+	    case XC_SAWTOOTH:
+		call strcpy ("sawtooth", XC_PSTRING(xc), SZ_FNAME)
+#	    case XC_MARK:
+#		call strcpy ("mark", XC_PSTRING(xc), SZ_FNAME)
+	    default:
+		;
+	    }
+	case NREFPTS:
+	    XC_NREFPTS(xc) = value
+	case CREGION:
+	    XC_CREGION(xc) = value
+	case NREGIONS:
+	    XC_NREGIONS(xc) = value
+	default:
+	    call error (0, "RG_XSETI: Undefined integer parameter.")
+	}
+end
+
+
+# RG_XSETP -- Set the value of a pointer parameter.
+
+procedure rg_xsetp (xc, param, value)
+
+pointer	xc		#I pointer to the cross-correlation structure
+int	param		#I parameter to be set
+pointer value		#O value of the pointer parameter
+
+begin
+	switch (param) {
+	case RC1:
+	    XC_RC1(xc) = value
+	case RC2:
+	    XC_RC2(xc) = value
+	case RL1:
+	    XC_RL1(xc) = value
+	case RL2:
+	    XC_RL2(xc) = value
+	case RZERO:
+	    XC_RZERO(xc) = value
+	case RXSLOPE:
+	    XC_RXSLOPE(xc) = value
+	case RYSLOPE:
+	    XC_RYSLOPE(xc) = value
+	case XSHIFTS:
+	    XC_XSHIFTS(xc) = value
+	case YSHIFTS:
+	    XC_YSHIFTS(xc) = value
+	case XCOR:
+	    XC_XCOR(xc) = value
+	case XREF:
+	    XC_XREF(xc) = value
+	case YREF:
+	    XC_YREF(xc) = value
+	case TRANSFORM:
+	    XC_TRANSFORM(xc) = value
+#	case CORAPODIZE:
+#	    XC_CORAPODIZE(xc) = value
+	default:
+	    call error (0, "RG_XSETP: Undefined pointer parameter.")
+	}
+end
+
+
+# RG_XSETR -- Set the value of a real parameter.
+
+procedure rg_xsetr (xc, param, value)
+
+pointer	xc		#I pointer to the cross-correlation structure
+int	param		#I parameter to be set
+real	value		#O value of real parameter
+
+begin
+	switch (param) {
+	case BVALUER:
+	    XC_BVALUER(xc) = value
+	case BVALUE:
+	    XC_BVALUE(xc) = value
+	case LOREJECT:
+	    XC_LOREJECT(xc) = value
+	case HIREJECT:
+	    XC_HIREJECT(xc) = value
+	case APODIZE:
+	    XC_APODIZE(xc) = value
+	case TXSHIFT:
+	    XC_TXSHIFT(xc) = value
+	case TYSHIFT:
+	    XC_TYSHIFT(xc) = value
+	default:
+	    call error (0, "RG_XSETR: Undefined real parameter.")
+	}
+end
+
+
+# RG_XSETS -- Set the value of a string parameter.
+
+procedure rg_xsets (xc, param, str)
+
+pointer	xc		#I pointer to the cross-correlation structure
+int	param		#I parameter to be set
+char	str[ARB]	#O value of string parameter
+
+int	index
+pointer	sp, temp
+int	strdic(), fnldir()
+
+begin
+	call smark (sp)
+	call salloc (temp, SZ_FNAME, TY_CHAR)
+
+	switch (param) {
+	case BSTRING:
+	    index = strdic (str, str, SZ_LINE, XC_BTYPES)
+	    if (index > 0) {
+	    	call strcpy (str, XC_BSTRING(xc), SZ_FNAME)
+		call rg_xseti (xc, BACKGRD, index)
+	    }
+	case FSTRING:
+	    index = strdic (str, str, SZ_LINE, XC_FTYPES)
+	    if (index > 0) {
+	        call strcpy (str, XC_FSTRING(xc), SZ_FNAME)
+		call rg_xseti (xc, FILTER, index)
+	    }
+	case CSTRING:
+	    index = strdic (str, str, SZ_LINE, XC_CTYPES)
+	    if (index > 0) {
+	        call strcpy (str, XC_CSTRING(xc), SZ_FNAME)
+		call rg_xseti (xc, CFUNC, index)
+	    }
+	case PSTRING:
+	    call strcpy (str, XC_PSTRING(xc), SZ_FNAME)
+	case REFIMAGE:
+	    call imgcluster (str, Memc[temp], SZ_FNAME)
+	    index = fnldir (Memc[temp], XC_REFIMAGE(xc), SZ_FNAME)
+	    call strcpy (Memc[temp+index], XC_REFIMAGE(xc), SZ_FNAME)
+	case IMAGE:
+	    call imgcluster (str, Memc[temp], SZ_FNAME)
+	    index = fnldir (Memc[temp], XC_IMAGE(xc), SZ_FNAME)
+	    call strcpy (Memc[temp+index], XC_IMAGE(xc), SZ_FNAME)
+	case OUTIMAGE:
+	    call strcpy (str, XC_OUTIMAGE(xc), SZ_FNAME)
+	case REGIONS:
+	    call strcpy (str, XC_REGIONS(xc), SZ_FNAME)
+	case DATABASE:
+	    index = fnldir (str, XC_DATABASE(xc), SZ_FNAME)
+	    call strcpy (str[index+1], XC_DATABASE(xc), SZ_FNAME)
+	case RECORD:
+	    call strcpy (str, XC_RECORD(xc), SZ_FNAME)
+	case REFFILE:
+	    index = fnldir (str, XC_REFFILE(xc), SZ_FNAME)
+	    call strcpy (str[index+1], XC_REFFILE(xc), SZ_FNAME)
+	default:
+	    call error (0, "RG_XSETS: Undefined string parameter.")
+	}
+
+	call sfree (sp)
+end
diff --git a/pkg/images/immatch/src/xregister/rgxtransform.x b/pkg/images/immatch/src/xregister/rgxtransform.x
new file mode 100644
index 00000000..63ee5f24
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxtransform.x
@@ -0,0 +1,446 @@
+include <imhdr.h>
+include <math.h>
+include "xregister.h"
+
+# RG_GXTRANSFORM -- Open the reference points file and the read the
+# coordinates of the reference points in the reference image. Return
+# the reference points file name and descriptor.
+
+int procedure rg_gxtransform (list, xc, reffile)
+
+int	list		#I list of reference points files
+pointer	xc		#I pointer to the cross-correlation structure
+char	reffile[ARB]	#O the output reference points file name
+
+int	tdf
+pointer	sp, line, pxref, pyref
+real	x1, y1, x2, y2, x3, y3
+int	fntgfnb(), open(), getline(), nscan()
+pointer	rg_xstatp()
+
+begin
+	# Get some working memory.
+	call smark (sp)
+	call salloc (line, SZ_LINE, TY_CHAR)
+
+	# Get the points to the reference point lists.
+	pxref = rg_xstatp (xc, XREF)
+	pyref = rg_xstatp (xc, YREF)
+	call aclrr (Memr[rg_xstatp(xc, XREF)], MAX_NREF)
+	call aclrr (Memr[rg_xstatp(xc, YREF)], MAX_NREF)
+
+	# Open the reference points file and read the coordinates.
+	while (fntgfnb (list, reffile, SZ_FNAME) != EOF) {
+
+	    iferr { 
+
+		# Open the reference file.
+	        tdf = open (reffile, READ_ONLY, TEXT_FILE) 
+		call aclrr (Memr[pxref], MAX_NREF)
+		call aclrr (Memr[pyref], MAX_NREF)
+
+		# Read up to three valid reference points from the list.
+	        while (getline (tdf, Memc[line]) != EOF) {
+		    call sscan (Memc[line])
+			call gargr (x1)
+			call gargr (y1)
+			call gargr (x2)
+			call gargr (y2)
+			call gargr (x3)
+			call gargr (y3)
+		    if (nscan () >= 2)
+			break
+		}
+
+		# Store the reference points.
+		if (nscan () == 2) {
+		    Memr[pxref] = x1
+		    Memr[pyref] = y1
+		    call rg_xseti (xc, NREFPTS, 1)
+		} else if (nscan () == 4) {
+		    Memr[pxref] = x1
+		    Memr[pyref] = y1
+		    Memr[pxref+1] = x2
+		    Memr[pyref+1] = y2
+		    call rg_xseti (xc, NREFPTS, 2)
+		} else if (nscan () == 6) {
+		    Memr[pxref] = x1
+		    Memr[pyref] = y1
+		    Memr[pxref+1] = x2
+		    Memr[pyref+1] = y2
+		    Memr[pxref+2] = x3
+		    Memr[pyref+2] = y3
+		    call rg_xseti (xc, NREFPTS, 2)
+		} else
+		    call rg_xseti (xc, NREFPTS, 0)
+
+	    } then { 
+		call rg_xseti (xc, NREFPTS, 0)
+	    }
+	}
+
+	call sfree (sp)
+
+	return (tdf)
+end
+
+
+# RG_ITRANSFORM -- Compute the transformation from the input image to the
+# reference image interactively.
+
+procedure rg_itransform (xc, imr, im, id)
+
+pointer	xc		#I pointer to the cross-correlation stucture
+pointer	imr		#I pointer to the reference image
+pointer	im		#I pointer to the input image
+pointer	id		#I pointer to the display device
+
+int	nref, nstar, wcs, key
+pointer	sp, cmd, x, y, pxref, pyref, ptrans
+real	wx, wy
+int	clgcur()
+pointer	rg_xstatp()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+	call salloc (x, MAX_NREF, TY_REAL)
+	call salloc (y, MAX_NREF, TY_REAL)
+	call aclrr (Memr[x], MAX_NREF)
+	call aclrr (Memr[y], MAX_NREF)
+
+	# Get the pointers.
+	pxref = rg_xstatp (xc, XREF)
+	pyref = rg_xstatp (xc, YREF)
+	ptrans = rg_xstatp (xc, TRANSFORM)
+
+	# Mark up to three reference stars.
+	nref = 0
+	call printf ("Mark reference star %d with the image cursor [q=quit]: ")
+	    call pargi (nref + 1)
+	while ((nref < MAX_NREF) && clgcur ("icommands", wx, wy, wcs, key,
+	    Memc[cmd], SZ_LINE) != EOF) {
+	    if (key == 'q') {
+		call printf ("\n")
+		break
+	    }
+	    if (wx < 0.5 || wx > IM_LEN(imr,1) + 0.5) {
+		call printf ("\n")
+		next
+	    }
+	    if (wy < 0.5 || wy > IM_LEN(imr,2) + 0.5) {
+		call printf ("\n")
+		next
+	    }
+	    call printf ("%g  %g\n")
+		call pargr (wx)
+		call pargr (wy)
+	    Memr[pxref+nref] = wx
+	    Memr[pyref+nref] = wy
+	    nref = nref + 1
+	    call rg_xseti (xc, NREFPTS, nref)
+	    if (nref >= MAX_NREF)
+		break
+	    call printf (
+	        "Mark reference star %d with the image cursor [q=quit]: ")
+	        call pargi (nref + 1)
+	}
+
+	# Mark the corresponding input image stars.
+	if (nref > 0) {
+
+	    nstar = 0
+	    call printf ("Mark image star %d with the image cursor [q=quit]: ")
+	        call pargi (nstar + 1)
+	    while ((nstar < nref) && clgcur ("icommands", wx, wy, wcs, key,
+	        Memc[cmd], SZ_LINE) != EOF) {
+	        if (key == 'q') {
+		    call printf ("\n")
+		    break
+		}
+	        if (wx < 0.5 || wx > IM_LEN(im,1) + 0.5) {
+		    call printf ("\n")
+		    next
+	        }
+	        if (wy < 0.5 || wy > IM_LEN(im,2) + 0.5) {
+		    call printf ("\n")
+		    next
+	        }
+	        call printf ("%g  %g\n")
+		    call pargr (wx)
+		    call pargr (wy)
+	        Memr[x+nstar] = wx
+	        Memr[y+nstar] = wy
+	        nstar = nstar + 1
+	        if (nstar >= MAX_NREF)
+		    break
+	        call printf (
+	            "Mark image star %d with the image cursor [q=quit]: ")
+	            call pargi (nstar + 1)
+	    }
+
+	    # Compute the transformation.
+	    if (nstar > 0) {
+	        switch (nstar) {
+	        case 0:
+	            call rg_xshift (Memr[pxref], Memr[pyref], Memr[pxref],
+	                Memr[pyref], Memr[ptrans])
+	        case 1:
+	            call rg_xshift (Memr[x], Memr[y], Memr[pxref], Memr[pyref], 
+	                Memr[ptrans])
+	        #case 2:
+	            #call rg_xtwostar (Memr[x], Memr[y], Memr[pxref],
+		        #Memr[pyref], Memr[ptrans])
+	        #case 3:
+	            #call rg_xthreestar (Memr[x], Memr[y], Memr[pxref],
+		        #Memr[pyref], Memr[ptrans])
+
+		default:
+	            call rg_xshift (Memr[pxref], Memr[pyref], Memr[pxref],
+	                Memr[pyref], Memr[ptrans])
+	        }
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# RG_XTRANSFORM -- Compute the transformation from the input image to
+# the reference image
+
+procedure rg_xtransform (tfd, xc)
+
+int	tfd		#I the reference points file descriptor
+pointer	xc		#I the cross-correlation file descriptor
+
+int	nref
+pointer	sp, line, x, y, pxref, pyref, ptrans
+int	getline(), rg_xstati(), nscan()
+pointer	rg_xstatp()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (line, SZ_LINE, TY_CHAR)
+	call salloc (x, MAX_NREF, TY_REAL)
+	call salloc (y, MAX_NREF, TY_REAL)
+	call aclrr (Memr[x], MAX_NREF)
+	call aclrr (Memr[y], MAX_NREF)
+
+	# Get the pointers to the reference image data.
+	nref = rg_xstati (xc, NREFPTS)
+	pxref = rg_xstatp (xc, XREF)
+	pyref = rg_xstatp (xc, YREF)
+	ptrans = rg_xstatp (xc, TRANSFORM)
+
+	# Read the input image reference points.
+	while ((nref > 0) && getline (tfd, Memc[line]) != EOF) {
+	    call sscan (Memc[line])
+		call gargr (Memr[x])
+		call gargr (Memr[y])
+		call gargr (Memr[x+1])
+		call gargr (Memr[y+1])
+		call gargr (Memr[x+2])
+		call gargr (Memr[y+2])
+	    if (nscan() >= 2 * nref)
+		break
+	}
+
+	# Compute the transform.
+	if (nscan () < 2 * nref) {
+	    call rg_xshift (Memr[pxref], Memr[pyref], Memr[pxref], Memr[pyref],
+	        Memr[ptrans])
+	} else {
+	    switch (nref) {
+	    case 0:
+	        call rg_xshift (Memr[pxref], Memr[pyref], Memr[pxref],
+	            Memr[pyref], Memr[ptrans])
+	    case 1:
+	        call rg_xshift (Memr[x], Memr[y], Memr[pxref], Memr[pyref], 
+	            Memr[ptrans])
+	    case 2:
+	        call rg_xtwostar (Memr[x], Memr[y], Memr[pxref], Memr[pyref],
+		    Memr[ptrans])
+	    case 3:
+	        call rg_xthreestar (Memr[x], Memr[y], Memr[pxref], Memr[pyref],
+		    Memr[ptrans])
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# RG_ETRANSFORM -- Evaulate the current transform at a single point.
+
+procedure rg_etransform (xc, xin, yin, xout, yout)
+
+pointer	xc		#I pointer to the cross-correlation structure
+real	xin, yin	#I the input x and y values
+real	xout, yout	#O the output x and y values
+
+pointer	ptrans
+pointer	rg_xstatp
+
+begin
+	ptrans = rg_xstatp (xc, TRANSFORM)
+	xout = Memr[ptrans] * xin + Memr[ptrans+1] * yin + Memr[ptrans+2]
+	yout = Memr[ptrans+3] * xin + Memr[ptrans+4] * yin + Memr[ptrans+5]
+end
+
+
+# RG_XSHIFT -- Compute the transformation coefficients required to define a
+# simple shift using a single data point.
+
+procedure rg_xshift (xref, yref, xlist, ylist, coeff)
+
+real	xref[ARB]		#I x reference coordinates
+real	yref[ARB]		#I y reference coordinates
+real	xlist[ARB]		#I x input coordinates
+real	ylist[ARB]		#I y input coordinates
+real	coeff[ARB]		#O output coefficient array
+
+begin
+	# Compute the x transformation.
+	coeff[1] = 1.0
+	coeff[2] = 0.0
+	coeff[3] = xref[1] - xlist[1]
+
+	# Compute the y transformation.
+	coeff[4] = 0.0
+	coeff[5] = 1.0
+	coeff[6] = yref[1] - ylist[1]
+end
+
+
+# RG_XTWOSTAR -- Compute the transformation coefficients required to
+# define a simple shift, magnification which is the same in x and y,
+# and rotation using two data points.
+
+procedure rg_xtwostar (xref, yref, xlist, ylist, coeff)
+
+real	xref[ARB]		#I x reference coordinates
+real	yref[ARB]		#I y reference coordinates
+real	xlist[ARB]		#I x input coordinates
+real	ylist[ARB]		#I y input coordinates
+real	coeff[ARB]		#O coefficient array
+
+real	rot, mag, dxlis, dylis, dxref, dyref, cosrot, sinrot
+real	rg_xposangle()
+
+begin
+	# Compute the deltas.
+	dxlis = xlist[2] - xlist[1]
+	dylis = ylist[2] - ylist[1]
+	dxref = xref[2] - xref[1]
+	dyref = yref[2] - yref[1]
+
+	# Compute the required rotation angle.
+	rot = rg_xposangle (dxref, dyref) - rg_xposangle (dxlis, dylis)
+	cosrot = cos (rot)
+	sinrot = sin (rot)
+
+	# Compute the required magnification factor.
+	mag = dxlis ** 2 + dylis ** 2
+	if (mag <= 0.0)
+	    mag = 0.0
+	else
+	    mag = sqrt ((dxref ** 2 + dyref ** 2) / mag)
+
+	# Compute the transformation coefficicents.
+	coeff[1] = mag * cosrot
+	coeff[2] = - mag * sinrot
+	coeff[3] = xref[1] - mag * cosrot * xlist[1] + mag * sinrot * ylist[1]
+	coeff[4] = mag * sinrot
+	coeff[5] = mag * cosrot
+	coeff[6] = yref[1] - mag * sinrot * xlist[1] - mag * cosrot * ylist[1]
+end
+
+
+# RG_THREESTAR -- Compute the transformation coefficients required to define
+# x and y shifts, x and ymagnifications, a rotation and skew, and a possible
+# axis flip using three tie points.
+
+procedure rg_xthreestar (xref, yref, xlist, ylist, coeff)
+
+real	xref[ARB]		#I x reference coordinates
+real	yref[ARB]		#I y reference coordinates
+real	xlist[ARB]		#I x input coordinates
+real	ylist[ARB]		#I y input coordinates
+real	coeff[ARB]		#O coefficient array
+
+real	dx23, dx13, dx12, dy23, dy13, dy12, det
+bool	fp_equalr()
+
+begin
+	# Compute the deltas.
+	dx23 = xlist[2] - xlist[3]
+	dx13 = xlist[1] - xlist[3]
+	dx12 = xlist[1] - xlist[2]
+	dy23 = ylist[2] - ylist[3]
+	dy13 = ylist[1] - ylist[3]
+	dy12 = ylist[1] - ylist[2]
+
+	# Compute the determinant.
+	det = xlist[1] * dy23 - xlist[2] * dy13 + xlist[3] * dy12
+	if (fp_equalr (det, 0.0)) {
+	    call rg_xtwostar (xref, yref, xlist, ylist, coeff)
+	    return
+	}
+
+	# Compute the x transformation.
+	coeff[1] = (xref[1] * dy23 - xref[2] * dy13 + xref[3] * dy12) / det
+	coeff[2] = (-xref[1] * dx23 + xref[2] * dx13 - xref[3] * dx12) / det
+	coeff[3] = (xref[1] * (xlist[2] * ylist[3] - xlist[3] * ylist[2]) +
+	    xref[2] * (ylist[1] * xlist[3] - xlist[1] * ylist[3]) +
+	    xref[3] * (xlist[1] * ylist[2] - ylist[1] * xlist[2])) / det
+
+	# Compute the y transformation.
+	coeff[4] = (yref[1] * dy23 - yref[2] * dy13 + yref[3] * dy12) / det
+	coeff[5] = (-yref[1] * dx23 + yref[2] * dx13 - yref[3] * dx12) / det
+	coeff[6] = (yref[1] * (xlist[2] * ylist[3] - xlist[3] * ylist[2]) +
+	    yref[2] * (ylist[1] * xlist[3] - xlist[1] * ylist[3]) +
+	    yref[3] * (xlist[1] * ylist[2] - ylist[1] * xlist[2])) / det
+end
+
+
+# RG_XPOSANGLE -- Compute the position angle of a 2D vector. The angle is
+# measured counter-clockwise from the positive x axis.
+
+real procedure rg_xposangle (x, y)
+
+real	x		#I x vector component
+real	y		#I y vector component
+
+real	theta
+bool	fp_equalr()
+
+begin
+	if (fp_equalr (y, 0.0)) {
+	    if (x > 0.0)
+		theta = 0.0
+	    else if (x < 0.0)
+		theta = PI
+	    else
+		theta = 0.0
+	} else if (fp_equalr (x, 0.0)) {
+	    if (y > 0.0)
+		theta = PI / 2.0
+	    else if (y < 0.0)
+		theta = 3.0 * PI / 2.0
+	    else
+		theta = 0.0
+	} else if (x > 0.0 && y > 0.0) {	# 1st quadrant
+	    theta = atan (y / x)
+	} else if (x > 0.0 && y < 0.0) {	# 4th quadrant
+	    theta = 2.0 * PI + atan (y / x)
+	} else if (x < 0.0 && y > 0.0) {	# 2nd quadrant
+	    theta = PI + atan (y / x)
+	} else if (x < 0.0 && y < 0.0) {	# 3rd quadrant
+	    theta = PI + atan (y / x)
+	}
+
+	return (theta)
+end
diff --git a/pkg/images/immatch/src/xregister/t_xregister.x b/pkg/images/immatch/src/xregister/t_xregister.x
new file mode 100644
index 00000000..f9fc9b22
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/t_xregister.x
@@ -0,0 +1,440 @@
+include <imhdr.h>
+include <fset.h>
+include <gset.h>
+include <imset.h>
+include "xregister.h"
+
+# T_XREGISTER -- Register a list of images using cross-correlation techniques.
+
+procedure t_xregister()
+
+pointer	freglist		# reference regions list
+pointer	database		# the shifts database
+int	dformat			# use the database format for the shifts file ?
+int	interactive		# interactive mode ?
+int	verbose			# verbose mode
+pointer	interpstr		# interpolant type
+int	boundary		# boundary extension type
+real	constant		# constant for boundary extension
+
+int	list1, listr, list2, reglist, reflist, reclist, tfd, stat, nregions
+int	c1, c2, l1, l2, ncols, nlines
+pointer	sp, image1, image2, imtemp, str, coords
+pointer	gd, id, imr, im1, im2, sdb, xc, mw
+real	shifts[2]
+bool	clgetb()
+int	imtopen(), imtlen(), imtgetim(), fntopnb(), clgwrd(), btoi()
+int	rg_xregions(), fntlenb(), rg_gxtransform(), rg_xstati()
+int	rg_xcorr(), rg_xicorr(), fntgfnb(), access(), open()
+pointer	gopen(), immap(), dtmap(), mw_openim()
+real	clgetr(), rg_xstatr()
+errchk	fntopnb(), gopen()
+
+begin
+	# Set STDOUT to flush on a newline character
+	call fseti (STDOUT, F_FLUSHNL, YES)
+
+	# Allocate temporary working space.
+	call smark (sp)
+
+	call salloc (freglist, SZ_LINE, TY_CHAR)
+	call salloc (image1, SZ_FNAME, TY_CHAR)
+	call salloc (image2, SZ_FNAME, TY_CHAR)
+	call salloc (imtemp, SZ_FNAME, TY_CHAR)
+	call salloc (database, SZ_FNAME, TY_CHAR)
+	call salloc (coords, SZ_FNAME, TY_CHAR)
+	call salloc (interpstr, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Get task parameters and open lists.
+	call clgstr ("input", Memc[str], SZ_LINE)
+	list1 = imtopen (Memc[str])
+	call clgstr ("reference", Memc[str], SZ_LINE)
+	listr = imtopen (Memc[str])
+	call clgstr ("regions", Memc[freglist], SZ_LINE)
+	call clgstr ("shifts", Memc[database], SZ_FNAME)
+	call clgstr ("output", Memc[str], SZ_LINE)
+	list2 = imtopen (Memc[str])
+	call clgstr ("records", Memc[str], SZ_LINE)
+	if (Memc[str] == EOS)
+	    reclist = NULL
+	else
+	    reclist = fntopnb (Memc[str], NO)
+	call clgstr ("coords", Memc[coords], SZ_LINE)
+
+	# Open the cross correlation fitting structure.
+	call rg_xgpars (xc)
+
+	# Test the reference image list length.
+	if (rg_xstati (xc, CFUNC) != XC_FILE) {
+	    if (imtlen (listr) <= 0)
+	        call error (0, "The reference image list is empty.")
+	    if (imtlen (listr) > 1 && imtlen (listr) != imtlen (list1))
+	        call error (0,
+	            "The number of reference and input images is not the same.")
+	    if (Memc[coords] == EOS)
+	        reflist = NULL
+	    else {
+	        reflist = fntopnb (Memc[coords], NO)
+		if (imtlen (listr) != fntlenb (reflist))
+	            call error (0,
+	      "The number of reference point files and images is not the same.")
+	   }
+	    iferr {
+	        reglist = fntopnb (Memc[freglist], NO)
+	    } then {
+	        reglist = NULL
+	    }
+	    call rg_xsets (xc, REGIONS, Memc[freglist])
+
+	} else {
+	    call imtclose (listr)
+	    listr = NULL
+	    reflist = NULL
+	    reglist = NULL
+	    call rg_xsets (xc, REGIONS, "")
+	}
+
+	# Close the output image list if it is empty.
+	if (imtlen (list2) == 0) {
+	    call imtclose (list2)
+	    list2 = NULL
+	}
+
+	# Check that the output image list is the same size as the input
+	# image list.
+	if (list2 != NULL) {
+	    if (imtlen (list1) != imtlen (list2)) {
+	        call imtclose (list1)
+	        if (list2 != NULL)
+	            call imtclose (list2)
+	       call error (0,
+	           "The number of input and output images is not the same.")
+	    }
+	}
+
+	# Check that the record list is the same length as the input
+	# image list length.
+	if (reclist != NULL) {
+	    if (fntlenb (reclist) != imtlen (list1))
+	        call error (0,
+	            "Input image and record lists are not the same length.")
+	}
+
+
+	# Open the database file.
+	dformat = btoi (clgetb ("databasefmt"))
+	if (rg_xstati (xc, CFUNC) == XC_FILE) {
+	    if (dformat == YES)
+	        sdb = dtmap (Memc[database], READ_ONLY)
+	    else
+		sdb = open (Memc[database], READ_ONLY, TEXT_FILE)
+	} else if (clgetb ("append")) {
+	    if (dformat == YES)
+	        sdb = dtmap (Memc[database], APPEND)
+	    else
+		sdb = open (Memc[database], NEW_FILE, TEXT_FILE)
+	} else if (access (Memc[database], 0, 0) == YES) {
+	    call error (0, "The shifts database file already exists")
+	} else {
+	    if (dformat == YES)
+	        sdb = dtmap (Memc[database], NEW_FILE)
+	    else
+		sdb = open (Memc[database], NEW_FILE, TEXT_FILE)
+	}
+	call rg_xsets (xc, DATABASE, Memc[database])
+
+	# Get the boundary extension parameters for the image shift.
+	call clgstr ("interp_type", Memc[interpstr], SZ_FNAME)
+	boundary = clgwrd ("boundary_type", Memc[str], SZ_LINE,
+	   "|constant|nearest|reflect|wrap|")
+	constant = clgetr ("constant")
+
+	if (rg_xstati (xc, CFUNC) == XC_FILE)
+	    interactive = NO
+	else
+	    interactive = btoi (clgetb ("interactive"))
+	if (interactive == YES) {
+	    call clgstr ("graphics", Memc[str], SZ_FNAME)
+	    iferr (gd = gopen (Memc[str], NEW_FILE, STDGRAPH))
+		gd = NULL
+	    call clgstr ("display", Memc[str], SZ_FNAME)
+	    iferr (id = gopen (Memc[str], APPEND, STDIMAGE))
+		id = NULL
+	    verbose = YES
+	} else {
+	    if (rg_xstati (xc, PFUNC) == XC_MARK)
+		call rg_xseti (xc, PFUNC, XC_CENTROID)
+	    gd = NULL
+	    id = NULL
+	    verbose = btoi (clgetb ("verbose"))
+	}
+
+	# Initialize the reference image filter descriptors
+	imr = NULL
+	tfd = NULL
+
+	# Initialize the overlap section.
+	c1 = INDEFI
+	c2 = INDEFI
+	l1 = INDEFI
+	l2 = INDEFI
+	ncols = INDEFI
+	nlines = INDEFI
+
+	# Do each set of input, reference,  and output images.
+	while ((imtgetim (list1, Memc[image1], SZ_FNAME) != EOF)) {
+
+	    # Open the reference image, and associated regions and coordinates
+	    # files if the correlation function is not file.
+
+	    if (rg_xstati (xc, CFUNC) != XC_FILE) {
+	        if (imtgetim (listr, Memc[str], SZ_FNAME) != EOF) {
+		    if (imr != NULL)
+		        call imunmap (imr)
+		    imr = immap (Memc[str], READ_ONLY, 0)
+		    if (IM_NDIM(imr) > 2)
+		        call error (0, "Reference images must be 1D or 2D")
+		    call rg_xsets (xc, REFIMAGE, Memc[str])
+		    nregions = rg_xregions (reglist, imr, xc, 1)
+		    if (nregions <= 0 && interactive == NO)
+		        call error (0, "The regions list is empty.")
+		    if (reflist != NULL) {
+		        if (tfd != NULL)
+			    call close (tfd)
+		        tfd = rg_gxtransform (reflist, xc, Memc[str])
+		        call rg_xsets (xc, REFFILE, Memc[str])
+		    }
+	        }
+	    } else
+		call rg_xsets (xc, REFIMAGE, "reference")
+
+	    # Open the input image.
+	    im1 = immap (Memc[image1], READ_ONLY, 0)
+	    if (IM_NDIM(im1) > 2) {
+	         call error (0, "Input images must be 1D or 2D")
+	    } else if (imr != NULL) {
+	        if (IM_NDIM(im1) != IM_NDIM(imr))
+		    call error (0,
+	      "Input images must have same dimensionality as reference images")
+	    }
+	    call imseti (im1, IM_TYBNDRY, BT_NEAREST)
+	    if (IM_NDIM(im1) == 1)
+	        call imseti (im1, IM_NBNDRYPIX, IM_LEN(im1,1))
+	    else
+	        call imseti (im1, IM_NBNDRYPIX,
+		    max (IM_LEN(im1,1), IM_LEN(im1,2)))
+	    call rg_xsets (xc, IMAGE, Memc[image1])
+
+	    # Open the output image if any.
+	    if (list2 == NULL) {
+	        im2 = NULL
+		Memc[image2] = EOS
+	    } else if (imtgetim (list2, Memc[image2], SZ_FNAME) != EOF) {
+		call xt_mkimtemp (Memc[image1], Memc[image2], Memc[imtemp],
+		    SZ_FNAME)
+	        im2 = immap (Memc[image2], NEW_COPY, im1)
+	    } else {
+	        im2 = NULL
+		Memc[image2] = EOS
+	    }
+	    call rg_xsets (xc, OUTIMAGE, Memc[image2])
+
+	    # Get the image record name for the shifts database.
+	    if (reclist == NULL)
+		call strcpy (Memc[image1], Memc[str], SZ_FNAME)
+	    else if (fntgfnb (reclist, Memc[str], SZ_FNAME) == EOF)
+		call strcpy (Memc[image1], Memc[str], SZ_FNAME)
+	    call rg_xsets (xc, RECORD, Memc[str])
+
+	    # Compute the initial coordinate shift.
+	    if (tfd != NULL)
+		call rg_xtransform (tfd, xc)
+
+	    # Perform the cross correlation function.
+	    if (interactive == YES) {
+		stat = rg_xicorr (imr, im1, im2, sdb, dformat, reglist, tfd,
+		    xc, gd, id)
+	    } else {
+	        stat = rg_xcorr (imr, im1, sdb, dformat, xc)
+		if (verbose == YES) {
+		    call rg_xstats (xc, REFIMAGE, Memc[str], SZ_LINE)
+	            call printf (
+		        "Average shift from %s to %s is %g %g pixels\n")
+	                call pargstr (Memc[image1])
+	                call pargstr (Memc[str])
+	                call pargr (rg_xstatr (xc, TXSHIFT))
+	                call pargr (rg_xstatr (xc, TYSHIFT))
+		}
+	    }
+
+	    # Compute the overlap region for the images.
+	    call rg_overlap (im1, rg_xstatr (xc, TXSHIFT),
+	        rg_xstatr (xc,TYSHIFT), c1, c2, l1, l2, ncols, nlines)
+
+	    # Shift the image and update the wcs.
+	    if (im2 != NULL && stat == NO) {
+		if (verbose == YES) {
+		    call printf (
+		        "\tShifting image %s to image %s ...\n")
+		        call pargstr (Memc[image1])
+			call pargstr (Memc[imtemp])
+		}
+
+		call rg_xshiftim (im1, im2, rg_xstatr (xc, TXSHIFT),
+		    rg_xstatr (xc, TYSHIFT), Memc[interpstr], boundary,
+		        constant)
+		mw = mw_openim (im1)
+		shifts[1] = rg_xstatr (xc, TXSHIFT)
+		shifts[2] = rg_xstatr (xc, TYSHIFT)
+		call mw_shift (mw, shifts, 03B)
+		call mw_saveim (mw, im2)
+		call mw_close (mw)
+	    }
+
+	    # Close up the input and output images.
+	    call imunmap (im1)
+	    if (im2 != NULL) {
+		call imunmap (im2)
+	        if (stat == YES)
+		    call imdelete (Memc[image2])
+		else
+		    call xt_delimtemp (Memc[image2], Memc[imtemp])
+	    }
+
+	    if (stat == YES)
+		break
+	    call rg_xindefr (xc)
+	}
+
+	if (verbose == YES)
+	    call rg_poverlap (c1, c2, l1, l2, ncols, nlines)
+
+	call rg_xfree (xc)
+
+	# Close up the lists.
+	if (imr != NULL)
+	    call imunmap (imr)
+	call imtclose (list1)
+	if (listr != NULL)
+	    call imtclose (listr)
+	if (reglist != NULL)
+	    call fntclsb (reglist)
+	if (list2 != NULL)
+	    call imtclose (list2)
+	if (tfd != NULL)
+	    call close (tfd)
+	if (reflist != NULL)
+	    call fntclsb (reflist)
+	if (reclist != NULL)
+	    call fntclsb (reclist)
+	if (dformat == YES)
+	    call dtunmap (sdb)
+	else
+	    call close (sdb)
+
+	# Close up the graphics and display devices.
+	if (gd != NULL)
+	    call gclose (gd)
+	if (id != NULL)
+	    call gclose (id)
+
+	call sfree (sp)
+end
+
+
+# RG_OVERLAP -- Compute the overlap region of the list of images.
+
+procedure rg_overlap (im1, xshift, yshift, x1, x2, y1, y2, ncols, nlines)
+
+pointer	im1			# pointer to the input image
+real	xshift			# the computed x shift of the input image
+real	yshift			# the computed y shift of the input image
+int	x1, x2			# the input/output column limits
+int	y1, y2			# the input/output line limits
+int	ncols, nlines		# the input/output size limits
+
+int	ixlo, ixhi, iylo, iyhi
+real	xlo, xhi, ylo, yhi
+
+begin
+	if (IS_INDEFR(xshift) || IS_INDEFR(yshift))
+	    return
+
+	# Compute the limits of the shifted image.
+	xlo = 1.0 + xshift
+	xhi = IM_LEN(im1,1) + xshift
+	ylo = 1.0 + yshift
+	yhi = IM_LEN(im1,2) + yshift
+
+	# Round up or down as appropriate.
+	ixlo = int (xlo)
+	if (xlo > ixlo)
+	    ixlo = ixlo + 1
+	ixhi = int (xhi)
+	if (xhi < ixhi)
+	    ixhi = ixhi - 1
+	iylo = int (ylo)
+	if (ylo > iylo)
+	    iylo = iylo + 1
+	iyhi = int (yhi)
+	if (yhi < iyhi)
+	    iyhi = iyhi - 1
+
+	# Determine the new limits.
+	if (IS_INDEFI(x1))
+	    x1 = ixlo
+	else
+	    x1 = max (ixlo, x1)
+	if (IS_INDEFI(x2))
+	    x2 = ixhi
+	else
+	    x2 = min (ixhi, x2)
+	if (IS_INDEFI(y1))
+	    y1 = iylo
+	else
+	    y1 = max (iylo, y1)
+	if (IS_INDEFI(y2))
+	    y2 = iyhi
+	else
+	    y2 = min (iyhi, y2)
+	if (IS_INDEFI(ncols))
+	    ncols = IM_LEN(im1,1)
+	else
+	    ncols = min (ncols, IM_LEN(im1,1))
+	if (IS_INDEFI(nlines))
+	    nlines = IM_LEN(im1,2)
+	else
+	    nlines = min (nlines, IM_LEN(im1,2))
+end
+
+
+# RG_POVERLAP -- Procedure to print the overlap and/or vignetted region.
+
+procedure rg_poverlap (x1, x2, y1, y2, ncols, nlines)
+
+int	x1, x2			# the input column limits
+int	y1, y2			# the input line limits
+int	ncols, nlines		# the number of lines and columns
+
+int	vx1, vx2, vy1, vy2
+
+begin
+	vx1 = max (1, min (x1, ncols))
+	vx2 = max (1, min (x2, ncols))
+	vy1 = max (1, min (y1, nlines))
+	vy2 = max (1, min (y2, nlines))
+
+	call printf ("Overlap region: [%d:%d,%d:%d]\n")
+	    call pargi (x1)
+	    call pargi (x2)
+	    call pargi (y1)
+	    call pargi (y2)
+	if (vx1 != x1 || vx2 != x2 || vy1 != y1 || vy2 != y2) {
+	    call printf ("Vignetted overlap region: [%d:%d,%d:%d]\n")
+	        call pargi (vx1)
+	        call pargi (vx2)
+	        call pargi (vy1)
+	        call pargi (vy2)
+	}
+end
diff --git a/pkg/images/immatch/src/xregister/xregister.h b/pkg/images/immatch/src/xregister/xregister.h
new file mode 100644
index 00000000..16c88b1e
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/xregister.h
@@ -0,0 +1,250 @@
+# Header file for XREGISTER
+
+# Define the cross correlation structure
+
+define	LEN_XCSTRUCT	(50 + 12 * SZ_FNAME + 12)
+
+define	XC_RC1		Memi[$1]	 # pointers to 1st column of ref regions
+define	XC_RC2		Memi[$1+1]	 # pointers to 2nd column of ref regions
+define	XC_RL1		Memi[$1+2]	 # pointers to 1st line of ref regions
+define	XC_RL2		Memi[$1+3]	 # pointers to 2nd line of ref regions
+define	XC_RZERO	Memi[$1+4]	 # pointers to zero pts of ref regions
+define	XC_RXSLOPE	Memi[$1+5]	 # pointers to x slopes of ref regions
+define	XC_RYSLOPE	Memi[$1+6]	 # pointers to y slopes of ref regions
+define	XC_XSHIFTS	Memi[$1+7]	 # pointers to x shifts of ref regions
+define	XC_YSHIFTS	Memi[$1+8]	 # pointers to y shifts of ref regions
+define	XC_NREGIONS	Memi[$1+9]	 # total number of regions
+define	XC_CREGION	Memi[$1+10]	 # the current region
+
+define	XC_NREFPTS	Memi[$1+11]	 # number of reference points
+define	XC_XREF		Memi[$1+12]	 # pointer to x reference points
+define	XC_YREF		Memi[$1+13]	 # pointer to y reference points
+define	XC_TRANSFORM	Memi[$1+14]	 # pointer to the transform
+define	XC_IXLAG	Memi[$1+15]	 # initial shift in x
+define	XC_IYLAG	Memi[$1+16]	 # initial shift in y
+define	XC_XLAG		Memi[$1+17]	 # current shift in x
+define	XC_YLAG		Memi[$1+18]	 # current shift in y
+define	XC_DXLAG	Memi[$1+19]	 # incremental shift in x
+define	XC_DYLAG	Memi[$1+20]	 # incremental shift in y
+
+define	XC_BACKGRD	Memi[$1+21]	 # type of background subtraction
+define	XC_BORDER	Memi[$1+22]	 # width of background border
+define	XC_BVALUER	Memr[P2R($1+23)] # reference background value
+define	XC_BVALUE	Memr[P2R($1+24)] # image bacground value
+define	XC_LOREJECT	Memr[P2R($1+25)] # low side rejection
+define	XC_HIREJECT	Memr[P2R($1+26)] # high side rejection
+define	XC_APODIZE	Memr[P2R($1+27)] # fraction of apodized region
+define	XC_FILTER	Memi[$1+28]	 # filter type
+
+define	XC_CFUNC	Memi[$1+30]	 # crosscor function
+define	XC_XWINDOW	Memi[$1+31]	 # width of correlation window in x
+define	XC_YWINDOW	Memi[$1+32]	 # width of correlation window in y
+define	XC_XCOR		Memi[$1+33]	 # pointer to cross-correlation function
+
+define	XC_PFUNC	Memi[$1+34]	 # correlation peak fitting function
+define	XC_XCBOX	Memi[$1+35]	 # x width of cor fitting box
+define	XC_YCBOX	Memi[$1+36]	 # y width of cor fitting box
+
+define	XC_TXSHIFT	Memr[P2R($1+37)] # total x shift
+define	XC_TYSHIFT	Memr[P2R($1+38)] # total y shift
+
+define	XC_BSTRING	Memc[P2C($1+50)]   	        # background type
+define	XC_FSTRING	Memc[P2C($1+50+SZ_FNAME+1)]     # filter string
+define	XC_CSTRING	Memc[P2C($1+50+2*SZ_FNAME+2)]   # cross-correlation type
+define	XC_PSTRING	Memc[P2C($1+50+3*SZ_FNAME+3)]   # peak centering
+
+define	XC_IMAGE	Memc[P2C($1+50+4*SZ_FNAME+4)]   # input image
+define	XC_REFIMAGE	Memc[P2C($1+50+5*SZ_FNAME+5)]   # reference image
+define	XC_REGIONS	Memc[P2C($1+50+6*SZ_FNAME+6)]   # regions list
+define	XC_DATABASE	Memc[P2C($1+50+7*SZ_FNAME+7)]   # shifts database
+define	XC_OUTIMAGE	Memc[P2C($1+50+8*SZ_FNAME+8)]   # output image
+define	XC_REFFILE	Memc[P2C($1+50+9*SZ_FNAME+9)]   # coordinates file
+define	XC_RECORD	Memc[P2C($1+50+10*SZ_FNAME+10)] # record
+
+# Define the id strings
+
+define	RC1		1
+define	RC2		2
+define	RL1		3
+define	RL2		4
+define	RZERO		5
+define	RXSLOPE		6
+define	RYSLOPE		7
+define	XSHIFTS		8
+define	YSHIFTS		9
+define	NREGIONS	10
+define	CREGION		11
+
+define	NREFPTS		12
+define	XREF		13
+define	YREF		14
+define	TRANSFORM	15
+define	IXLAG		16
+define	IYLAG		17
+define	XLAG		18
+define	YLAG		19
+define	DXLAG		20
+define	DYLAG		21
+
+define	BACKGRD		22
+define	BVALUER		23
+define	BVALUE		24
+define	BORDER		25
+define	LOREJECT	26
+define	HIREJECT	27
+define	APODIZE		28
+define	FILTER		29
+
+define	CFUNC		30
+define	XWINDOW		31
+define	YWINDOW		32
+define	XCOR		33
+
+define	PFUNC		34
+define	XCBOX		35
+define	YCBOX		36
+
+define	TXSHIFT		37
+define	TYSHIFT		38
+
+define	CSTRING		39
+define	BSTRING		40
+define	PSTRING		41
+define	FSTRING		42
+
+define	IMAGE		43
+define	REFIMAGE	44
+define	REGIONS		45
+define	OUTIMAGE	46
+define	REFFILE		47
+define	DATABASE	48
+define	RECORD		49
+
+# Define the default parameter values
+
+define	DEF_IXLAG	0
+define	DEF_IYLAG	0
+define	DEF_DXLAG	0
+define	DEF_DYLAG	0
+define	DEF_XWINDOW	5
+define	DEF_YWINDOW	5
+
+define	DEF_BACKGRD	XC_BNONE
+define	DEF_BORDER	INDEFI
+define  DEF_LOREJECT	INDEFR
+define	DEF_HIREJECT	INDEFR
+
+define	DEF_XCBOX	5
+define  DEF_YCBOX	5
+define	DEF_PFUNC	XC_CENTROID
+
+# Define the background fitting techniques
+
+define	XC_BNONE	1
+define	XC_MEAN		2
+define	XC_MEDIAN	3
+define	XC_SLOPE	4
+
+define  XC_BTYPES	"|none|mean|median|plane|"
+
+# Define the filtering options
+
+define	XC_FNONE	1
+define	XC_LAPLACE	2
+
+define  XC_FTYPES	"|none|laplace|"
+
+# Define the cross correlation techniques
+
+define	XC_DISCRETE	1
+define	XC_FOURIER	2
+define	XC_DIFFERENCE	3
+define	XC_FILE		4
+
+define	XC_CTYPES	"|discrete|fourier|difference|file|"
+
+# Define the peak fitting functions
+
+define	XC_PNONE	1
+define	XC_CENTROID	2
+define	XC_SAWTOOTH	3
+define	XC_PARABOLA	4
+define	XC_MARK		5
+
+define	XC_PTYPES	"|none|centroid|sawtooth|parabola|mark|"
+
+# Miscellaneous
+
+define	MAX_NREGIONS	100
+define	MAX_NREF	3
+define	MAX_NTRANSFORM	6
+
+# Commands
+
+define  XCMDS "|reference|input|regions|shifts|output|records|transform|\
+cregion|xlag|ylag|dxlag|dylag|background|border|loreject|hireject|apodize|\
+filter|correlation|xwindow|ywindow|function|xcbox|ycbox|show|mark|"
+
+define  XSHOW	"|data|background|correlation|center|"
+
+define	XSHOW_DATA			1
+define	XSHOW_BACKGROUND		2
+define	XSHOW_CORRELATION		3
+define	XSHOW_PEAKCENTER		4
+
+define	XCMD_REFIMAGE		1
+define	XCMD_IMAGE		2
+define	XCMD_REGIONS		3
+define	XCMD_DATABASE		4
+define	XCMD_OUTIMAGE		5
+define	XCMD_RECORD		6
+define	XCMD_REFFILE		7
+define	XCMD_CREGION		8
+define	XCMD_XLAG		9
+define	XCMD_YLAG		10
+define	XCMD_DXLAG		11
+define	XCMD_DYLAG		12
+define	XCMD_BACKGROUND		13
+define	XCMD_BORDER		14
+define	XCMD_LOREJECT		15
+define	XCMD_HIREJECT		16
+define	XCMD_APODIZE		17
+define	XCMD_FILTER		18
+define	XCMD_CORRELATION	19
+define	XCMD_XWINDOW		20
+define	XCMD_YWINDOW		21
+define	XCMD_PEAKCENTER		22
+define	XCMD_XCBOX		23
+define	XCMD_YCBOX		24
+define	XCMD_SHOW		25
+define	XCMD_MARK		26
+
+# Keywords
+
+define	KY_REFIMAGE		"reference"
+define	KY_IMAGE		"input"
+define	KY_REGIONS		"regions"
+define	KY_DATABASE		"shifts"
+define	KY_OUTIMAGE		"output"
+define	KY_RECORD		"record"
+define	KY_REFFILE		"coords"
+define	KY_NREGIONS		"nregions"
+define	KY_CREGION		"region"
+define	KY_XLAG			"xlag"
+define	KY_YLAG			"ylag"
+define	KY_DXLAG		"dxlag"
+define	KY_DYLAG		"dylag"
+define	KY_BACKGROUND		"background"
+define	KY_BORDER		"border"
+define	KY_LOREJECT		"loreject"
+define	KY_HIREJECT		"hireject"
+define	KY_APODIZE		"apodize"
+define	KY_FILTER		"filter"
+define	KY_CORRELATION		"correlation"
+define	KY_XWINDOW		"xwindow"
+define	KY_YWINDOW		"ywindow"
+define	KY_PEAKCENTER		"function"
+define	KY_XCBOX		"xcbox"
+define	KY_YCBOX		"ycbox"
+define	KY_TXSHIFT		"xshift"
+define	KY_TYSHIFT		"yshift"
diff --git a/pkg/images/immatch/src/xregister/xregister.key b/pkg/images/immatch/src/xregister/xregister.key
new file mode 100644
index 00000000..1956c88f
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/xregister.key
@@ -0,0 +1,47 @@
+		Interactive Keystroke Commands
+
+?	Print help 
+:	Colon commands
+t	Define the offset between the reference and input images
+c	Draw a contour plot of the cross-correlation function
+x	Draw a column plot of the cross-correlation function
+y	Draw a line plot of the cross-correlation function
+r	Redraw the current plot
+f	Recompute the cross-correlation function
+o	Enter the image overlay plot submenu 
+w	Update the task parameters
+q	Exit
+
+
+		Colon Commands
+
+:mark		    Mark regions on the display
+:show	            Show current values of all the parameters
+
+
+		Show/set Parameters
+
+:reference	[string]    Show/set the current reference image name
+:input		[string]    Show/set the current input image name
+:regions	[string]    Show/set the regions to be cross-correlated
+:shifts		{string]    Show/set the shifts database file name
+:coords		[string]    Show/set the current coordinates file name
+:output		[string]    Show/set the current output image name
+:records	[string]    Show/set the current database record name
+:xlag		[value]     Show/set the initial lag in x
+:ylag		[value]     Show/set the initial lag in y
+:dxlag		[value]     Show/set the incremental lag in x
+:dylag		[value]     Show/set the incremental lag in y
+:cregion	[value]	    Show/set the current region
+:background	[string]    Show/set the background fitting function
+:border		[value]     Show/set border region for background fitting
+:loreject	[value]     Show/set low side k-sigma rejection parameter
+:hireject	[value]     Show/set high side k-sigma rejection parameter
+:apodize	[value]	    Show/set percent of end points to apodize
+:filter		[string]    Show/set the default spatial filter 
+:correlation	[string]    Show/set the cross-correlation function 
+:xwindow	[value]     Show/set width of cross-correlation window in x
+:ywindow	[value]     Show/set width of cross-correlation window in y
+:function	[string]    Show/set correlation peak centering function 
+:xcbox		[value]	    Show/set the centering box width in x
+:ycbox		[value]	    Show/set the centering box width in y