Repatch (from linux) of OSX IRAF

1 files changed, 391 insertions, 0 deletions

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