Initial commit

1 files changed, 624 insertions, 0 deletions

diff --git a/pkg/images/imgeom/src/t_magnify.x b/pkg/images/imgeom/src/t_magnify.x
new file mode 100644
index 00000000..ed797500
--- /dev/null
+++ b/pkg/images/imgeom/src/t_magnify.x
@@ -0,0 +1,624 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <error.h>
+include	<imhdr.h>
+include	<imset.h>
+include	<math/iminterp.h>
+
+# T_MAGNIFY -- Change coordinate origin and pixel interval in 2D images.
+#
+# The input and output images are given by image template lists.  The
+# number of output images must match the number of input images.  Image
+# sections are allowed in the input images and ignored in the output
+# images.  If the input and output image names are the same then the
+# magnification is performed to a temporary file which then replaces
+# the input image.
+
+# Interpolation types and boundary extension types.
+
+define	BTYPES	"|constant|nearest|reflect|wrap|project|"
+define	SZ_BTYPE	8
+
+procedure t_magnify ()
+
+pointer	input				# Pointer to input image list
+pointer	output				# Pointer to output image list
+pointer	interp				# Pointer to image interpolation type
+pointer	boundary			# Pointer to boundary extension type
+real	bconst				# Boundary extension pixel value
+real	xmag, ymag			# Image magnifications
+real	dx, dy				# Step size
+real	x1, y1				# Starting coordinates
+real	x2, y2				# Ending coordinates
+int	flux				# Flux conserve
+
+int	list1, list2, btype, logfd
+pointer	sp, in, out, image1, image2, image3, mw, errmsg
+real	a, b, c, d, shifts[2], scale[2]
+
+bool	clgetb(), envgetb(), fp_equalr()
+int	clgwrd(), imtopen(), imtgetim(), imtlen(), open(), btoi(), errget()
+pointer	mw_openim(), immap()
+real	clgetr()
+errchk	open(), mg_magnify1(), mg_magnify2()
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_LINE, TY_CHAR)
+	call salloc (output, SZ_LINE, TY_CHAR)
+	call salloc (interp, SZ_FNAME, TY_CHAR)
+	call salloc (boundary, SZ_BTYPE, TY_CHAR)
+	call salloc (image1, SZ_FNAME, TY_CHAR)
+	call salloc (image2, SZ_FNAME, TY_CHAR)
+	call salloc (image3, SZ_FNAME, TY_CHAR)
+	call salloc (errmsg, SZ_FNAME, TY_CHAR)
+
+	# Get task parameters.
+	call clgstr ("input", Memc[input], SZ_LINE)
+	call clgstr ("output", Memc[output], SZ_LINE)
+	call clgstr ("interpolation", Memc[interp], SZ_FNAME)
+	btype = clgwrd ("boundary", Memc[boundary], SZ_BTYPE, BTYPES)
+	bconst = clgetr ("constant")
+	a = clgetr ("x1")
+	b = clgetr ("x2")
+	dx = clgetr ("dx")
+	c = clgetr ("y1")
+	d = clgetr ("y2")
+	dy = clgetr ("dy")
+	flux = btoi (clgetb ("fluxconserve")) 
+
+	# If the pixel interval INDEF then use the a magnification factor
+	# to determine the pixel interval.
+
+	if (IS_INDEF (dx)) {
+	    xmag = clgetr ("xmag")
+	    if (xmag < 0.0)
+		dx = -xmag
+	    else if (xmag > 0.0)
+		dx = 1.0 / xmag
+	    else
+		dx = 0.0
+	}
+
+	if (IS_INDEF (dy)) {
+	    ymag = clgetr ("ymag")
+	    if (ymag < 0.0)
+		dy = -ymag
+	    else if (ymag > 0.0)
+		dy = 1.0 / ymag
+	    else
+		dy = 0.0
+	}
+
+	if (fp_equalr (dx, 0.0) || fp_equalr (dy, 0.0)) {
+	    call error (0, "Illegal magnification")
+	} else {
+	    xmag = 1.0 / dx
+	    ymag = 1.0 / dy
+	}
+
+
+	# Open the log file.
+	call clgstr ("logfile", Memc[image1], SZ_FNAME)
+	iferr (logfd = open (Memc[image1], APPEND, TEXT_FILE))
+	    logfd = NULL
+
+	# Expand the input and output image lists.
+	list1 = imtopen (Memc[input])
+	list2 = imtopen (Memc[output])
+	if (imtlen (list1) != imtlen (list2)) {
+	    call imtclose (list1)
+	    call imtclose (list2)
+	    call error (0, "Number of input and output images not the same")
+	}
+
+	# Magnify each set of input/output images with the 2D interpolation
+	# package.
+
+	while ((imtgetim (list1, Memc[image1], SZ_FNAME) != EOF) &&
+	    (imtgetim (list2, Memc[image2], SZ_FNAME) != EOF)) {
+
+	    # Map the input and output images.
+	    call xt_mkimtemp (Memc[image1], Memc[image2], Memc[image3],
+	        SZ_FNAME)
+	    in = immap (Memc[image1], READ_ONLY, 0)
+	    out = immap (Memc[image2], NEW_COPY, in)
+
+	    # Set the limits of the output image.
+	    x1 = a
+	    x2 = b
+	    y1 = c
+	    y2 = d
+
+	    # Magnify the image making sure to update the wcs.
+	    iferr {
+	        if (IM_NDIM(in) == 1) {
+		    call mg_magnify1 (in, out, Memc[interp], btype, bconst,
+		        x1, x2, dx, flux)
+		    if (!envgetb ("nomwcs")) {
+		        mw = mw_openim (in)
+		        scale[1] = xmag
+		        shifts[1] =  1. - xmag * x1
+		        call mw_scale (mw, scale, 01B)
+		        call mw_shift (mw, shifts, 01B)
+		        call mw_saveim (mw, out)
+		        call mw_close (mw)
+		    }
+	        } else if (IM_NDIM(in) == 2) {
+	            call mg_magnify2 (in, out, Memc[interp], btype, bconst,
+		        x1, x2, dx, y1, y2, dy, flux)
+		    if (!envgetb ("nomwcs")) {
+		        mw = mw_openim (in)
+		        scale[1] = xmag
+		        scale[2] = ymag
+		        shifts[1] =  1. - xmag * x1
+		        shifts[2] =  1. - ymag * y1
+		        call mw_scale (mw, scale, 03B)
+		        call mw_shift (mw, shifts, 03B)
+		        call mw_saveim (mw, out)
+		        call mw_close (mw)
+		    }
+	        } else {
+	            call imunmap (out)
+	            call imunmap (in)
+	            call xt_delimtemp (Memc[image2], Memc[image3])
+		    if (logfd != NULL) {
+		        call fprintf (logfd, "\n%s\n")
+		            call pargstr (Memc[image3])
+		        call fprintf (logfd,
+		            "  Cannot magnify image %s to image %s.\n")
+		            call pargstr (Memc[image1])
+		            call pargstr (Memc[image3])
+		        call fprintf (logfd,
+			"  Dimensions are greater than 2.\n")
+		    }
+		    call imdelete (Memc[image3])
+		    next
+	        }
+
+	    } then {
+		 if (logfd != NULL) {
+		    call fprintf (logfd, "\n%s\n")
+		        call pargstr (Memc[image3])
+		    call fprintf (logfd,
+		        "  Cannot magnify image %s to image %s.\n")
+		         call pargstr (Memc[image1])
+		         call pargstr (Memc[image3])
+		    if (errget (Memc[errmsg], SZ_FNAME) <= 0)
+			;
+		    call fprintf (logfd, "%s")
+			call pargstr (Memc[errmsg])
+		 }
+		 call imdelete (Memc[image3])
+	         call imunmap (out)
+	         call imunmap (in)
+	         call xt_delimtemp (Memc[image2], Memc[image3])
+	    } else {
+
+	        if (logfd != NULL) {
+		    call fprintf (logfd, "\n%s\n")
+		        call pargstr (Memc[image3])
+		    call fprintf (logfd, "  Magnify image %s to image %s.\n")
+		        call pargstr (Memc[image1])
+		        call pargstr (Memc[image3])
+		    call fprintf (logfd, "  Interpolation is %s.\n")
+		        call pargstr (Memc[interp])
+		    call fprintf (logfd, "  Boundary extension is %s.\n")
+		        call pargstr (Memc[boundary])
+		    if (btype == 1) {
+		        call fprintf (logfd,
+			    "  Boundary pixel constant is %g.\n")
+			    call pargr (bconst)
+		    }
+		    call fprintf (logfd,
+		        "  Output coordinates in terms of input coordinates:\n")
+		    call fprintf (logfd,
+		        "    x1 = %10.4g, x2 = %10.4g, dx = %10.6g\n")
+		        call pargr (x1)
+		        call pargr (x2)
+		        call pargr (dx)
+		    if (IM_NDIM(in) == 2) {
+		        call fprintf (logfd,
+		            "    y1 = %10.4g, y2 = %10.4g, dy = %10.6g\n")
+		            call pargr (y1)
+		            call pargr (y2)
+		            call pargr (dy)
+		    }
+	        }
+
+	        call imunmap (out)
+	        call imunmap (in)
+	        call xt_delimtemp (Memc[image2], Memc[image3])
+	    }
+
+	}
+
+	call imtclose (list1)
+	call imtclose (list2)
+	call close (logfd)
+	call sfree (sp)
+end
+
+
+define	NYOUT2		16	# Number of input lines to use for interpolation
+define	NMARGIN		3	# Number of edge lines to add for interpolation
+define	NMARGIN_SPLINE3	16	# Number of edge lines to add for interpolation
+
+
+# MG_MAGNIFY1 -- Magnify the input input image to create the output image.
+
+procedure mg_magnify1 (in, out, interp, btype, bconst, x1, x2, dx, flux)
+
+pointer	in			# pointer to the input image
+pointer	out			# pointer to the output image
+char	interp[ARB]		# Interpolation type
+int	btype			# Boundary extension type
+real	bconst			# Boundary extension constant
+real	x1, x2			# Starting and ending points of output image
+real	dx			# Pixel interval
+int	flux			# Conserve flux?
+
+int	i, nxin, nxout, nxymargin, itype, nsinc, nincr, col1, col2
+pointer	sp, x, z, buf, asi
+real	xshift
+pointer	imgs1r(), impl1r()
+int	asigeti()
+
+begin
+	# Set the default values for the output image limits if they are INDEF
+	# and calculate the number of output pixels.
+
+	if (IS_INDEF (x1))
+	    x1 = 1.
+	if (IS_INDEF (x2))
+	    x2 = IM_LEN (in, 1)
+	if (x1 > x2)
+	    call error (0, "  X1 cannot be greater than X2\n")
+
+	# Set the number of output pixels in the image header.
+
+	nxout = (x2 - x1) / dx + 1
+	IM_LEN(out, 1) = nxout
+
+	# Initialize the interpolator.
+
+	call asitype (interp, itype, nsinc, nincr, xshift)
+	call asisinit (asi, itype, nsinc, nincr, xshift, 0.0)
+
+	# Round the coordinate limits to include the output image coordinate
+	# limits and the set boundary.
+
+	col1 = x1
+	col2 = nint (x2)
+	if (itype == II_SPLINE3)
+	    nxymargin = NMARGIN_SPLINE3
+	else if (itype == II_SINC || itype == II_LSINC)
+	    nxymargin = asigeti (asi, II_ASINSINC)
+	else if (itype == II_DRIZZLE)
+	    nxymargin = max (nint (dx), NMARGIN)
+	else
+	    nxymargin = NMARGIN
+	call mg_setboundary1 (in, col1, col2, btype, bconst, nxymargin)
+	col1 = col1 - nxymargin
+	if (col1 <= 0)
+	    col1 = col1 - 1
+	col2 = col2 + nxymargin + 1
+
+	# Allocate memory for the interpolation coordinates.
+	# Also initialize the image data buffer.
+
+	call smark (sp)
+	call salloc (x, 2 * nxout, TY_REAL)
+
+	# Set the x interpolation coordinates.  The coordinates are relative
+	# to the boundary extended input image.
+
+	if (itype == II_DRIZZLE) {
+	    do i = 1, nxout {
+	        Memr[x+2*i-2] = x1 + (i - 1.5) * dx - col1 + 1
+	        Memr[x+2*i-1] = x1 + (i - 0.5) * dx - col1 + 1
+	    }
+	} else {
+	    do i = 1, nxout
+	        Memr[x+i-1] = x1 + (i - 1) * dx - col1 + 1
+	}
+
+	# Fit the output image.
+	nxin = col2 - col1 + 1
+	buf = imgs1r (in, col1, col2)
+	call asifit (asi, Memr[buf], nxin)
+
+	# Evaluate the output image pixel values.
+	z = impl1r (out)
+	call asivector (asi, Memr[x], Memr[z], nxout)
+	#if (itype != II_DRIZZLE && flux == YES)
+	if (flux == YES)
+	    call amulkr (Memr[z], dx, Memr[z], nxout)
+
+	# Free memory and unmap the images.
+	call asifree (asi)
+	call sfree (sp)
+end
+
+
+# MG_MAGNIFY2 -- Magnify the input input image to create the output image.
+
+procedure mg_magnify2 (in, out, interp, btype, bconst, x1, x2, dx, y1, y2,
+        dy, flux)
+
+pointer	in			# pointer to the input image
+pointer	out			# pointer to the output image
+char	interp[ARB]		# Interpolation type
+int	btype			# Boundary extension type
+real	bconst			# Boundary extension constant
+real	x1, y1			# Starting point of output image
+real	x2, y2			# Ending point of output image
+real	dx, dy			# Pixel interval
+int	flux			# Conserve flux?
+
+int	i, nxin, nxout, nyout, nxymargin, itype, nsinc, nincr
+int	l1out, l2out, nlout, l1in, l2in, nlin, fstline, lstline
+int	col1, col2, line1, line2
+real	shift
+pointer	msi
+pointer	sp, x, y, z, buf
+
+pointer	imps2r()
+int	msigeti()
+
+begin
+	# Set the default values for the output image limits if they are INDEF
+	# and calculate the number of output pixels.
+
+	if (IS_INDEF (x1))
+	    x1 = 1.
+	if (IS_INDEF (x2))
+	    x2 = IM_LEN (in, 1)
+	if (IS_INDEF (y1))
+	    y1 = 1.
+	if (IS_INDEF (y2))
+	    y2 = IM_LEN (in, 2)
+	if (x1 > x2)
+	    call error (0, "  X1 cannot be greater than X2\n")
+	if (y1 > y2)
+	    call error (0, "  Y1 cannot be greater than Y2\n")
+	nxout = (x2 - x1) / dx + 1
+	nyout = (y2 - y1) / dy + 1
+
+	# Set the number of output pixels in the image header.
+
+	IM_LEN(out, 1) = nxout
+	IM_LEN(out, 2) = nyout
+
+	# Initialize the interpolator.
+
+	call msitype (interp, itype, nsinc, nincr, shift)
+	call msisinit (msi, itype, nsinc, nincr, nincr, shift, shift, 0.0)
+
+	# Compute the number of margin pixels required
+
+	if (itype == II_BISPLINE3)
+	    nxymargin = NMARGIN_SPLINE3
+	else if (itype == II_BISINC || itype == II_BILSINC)
+	    nxymargin = msigeti (msi, II_MSINSINC)
+	else if (itype == II_BIDRIZZLE)
+	    nxymargin = max (nint (dx), nint(dy), NMARGIN)
+	else
+	    nxymargin = NMARGIN
+
+	# Round the coordinate limits to include the output image coordinate
+	# limits and the set boundary.
+
+	col1 = x1
+	col2 = nint (x2)
+	line1 = y1
+	line2 = nint (y2)
+	call mg_setboundary2 (in, col1, col2, line1, line2, btype, bconst,
+	    nxymargin)
+
+	# Compute the input image column limits.
+	col1 = col1 - nxymargin
+	if (col1 <= 0)
+	    col1 = col1 - 1
+	col2 = col2 + nxymargin + 1
+	nxin = col2 - col1 + 1
+
+	# Allocate memory for the interpolation coordinates.
+	# Also initialize the image data buffer.
+
+	call smark (sp)
+	call salloc (x, 2 * nxout, TY_REAL)
+	call salloc (y, 2 * NYOUT2, TY_REAL)
+	buf = NULL
+	fstline = 0
+	lstline = 0
+
+	# Set the x interpolation coordinates which do not change from
+	# line to line.  The coordinates are relative to the boundary
+	# extended input image.
+
+	if (itype == II_BIDRIZZLE) {
+	    do i = 1, nxout {
+	        Memr[x+2*i-2] = x1 + (i - 1.5) * dx - col1 + 1
+	        Memr[x+2*i-1] = x1 + (i - 0.5) * dx - col1 + 1
+	    }
+	} else {
+	    do i = 1, nxout
+	        Memr[x+i-1] = x1 + (i - 1) * dx - col1 + 1
+	}
+
+	# Loop over the image sections.
+	for (l1out = 1; l1out <= nyout; l1out = l1out + NYOUT2) {
+
+	    # Define the range of output lines.
+	    l2out = min (l1out + NYOUT2 - 1, nyout)
+	    nlout = l2out - l1out + 1
+
+	    # Define the corresponding range of input lines.
+	    l1in = y1 + (l1out - 1) * dy - nxymargin
+	    if (l1in <= 0)
+		l1in = l1in - 1
+	    l2in = y1 + (l2out - 1) * dy + nxymargin + 1
+	    nlin = l2in - l1in + 1
+
+	    # Get the apporiate image section and compute the coefficients.
+	    if ((buf == NULL) || (l1in < fstline) || (l2in > lstline)) {
+		fstline = l1in
+		lstline = l2in
+		call mg_bufl2r (in, col1, col2, l1in, l2in, buf)
+		call msifit (msi, Memr[buf], nxin, nlin, nxin)
+	    }
+
+	    # Output the section.
+	    z = imps2r (out, 1, nxout, l1out, l2out)
+
+	    # Compute the y values.
+	    if (itype == II_BIDRIZZLE) {
+	        do i = l1out, l2out {
+		    Memr[y+2*(i-l1out)] = y1 + (i - 1.5) * dy - fstline + 1
+		    Memr[y+2*(i-l1out)+1] = y1 + (i - 0.5) * dy - fstline + 1
+		}
+	    } else {
+	        do i = l1out, l2out
+		    Memr[y+i-l1out] = y1 + (i - 1) * dy - fstline + 1
+	    }
+
+	    # Evaluate the interpolant.
+	    call msigrid (msi, Memr[x], Memr[y], Memr[z], nxout, nlout, nxout)
+	    if (flux == YES)
+		call amulkr (Memr[z], dx * dy, Memr[z], nxout * nlout)
+	}
+
+	# Free memory and buffers.
+	call msifree (msi)
+	call mfree (buf, TY_REAL)
+	call sfree (sp)
+end
+
+
+# MG_BUFL2R -- Maintain buffer of image lines.  A new buffer is created when
+# the buffer pointer is null or if the number of lines requested is changed.
+# The minimum number of image reads is used.
+
+procedure mg_bufl2r (im, col1, col2, line1, line2, buf)
+
+pointer	im		# Image pointer
+int	col1		# First image column of buffer
+int	col2		# Last image column of buffer
+int	line1		# First image line of buffer
+int	line2		# Last image line of buffer
+pointer	buf		# Buffer
+
+int	i, ncols, nlines, nclast, llast1, llast2, nllast
+pointer	buf1, buf2
+
+pointer	imgs2r()
+
+begin
+	ncols = col2 - col1 + 1
+	nlines = line2 - line1 + 1
+
+	# If the buffer pointer is undefined then allocate memory for the
+	# buffer.  If the number of columns or lines requested changes
+	# reallocate the buffer.  Initialize the last line values to force
+	# a full buffer image read.
+
+	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
+	}
+
+	# Read only the image lines with are different from the last buffer.
+
+	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)
+	    }
+	}
+
+	# Save the buffer parameters.
+
+	llast1 = line1
+	llast2 = line2
+	nclast = ncols
+	nllast = nlines
+end
+
+
+# MG_SETBOUNDARY1 -- Set boundary extension for a 1D image.
+
+procedure mg_setboundary1 (im, col1, col2, btype, bconst, nxymargin)
+
+pointer	im			# IMIO pointer
+int	col1, col2		# Range of columns
+int	btype			# Boundary extension type
+real	bconst			# Constant for constant boundary extension
+int	nxymargin		# Number of margin pixels
+
+int	btypes[5]
+int	nbndrypix
+
+data	btypes /BT_CONSTANT, BT_NEAREST, BT_REFLECT, BT_WRAP, BT_PROJECT/
+
+begin
+	nbndrypix = 0
+	nbndrypix = max (nbndrypix, 1 - col1)
+	nbndrypix = max (nbndrypix, col2 - IM_LEN(im, 1))
+
+	call imseti (im, IM_TYBNDRY, btypes[btype])
+	call imseti (im, IM_NBNDRYPIX, nbndrypix + nxymargin + 1)
+	if (btypes[btype] == BT_CONSTANT)
+	    call imsetr (im, IM_BNDRYPIXVAL, bconst)
+end
+
+
+# MG_SETBOUNDARY2 -- Set boundary extension for a 2D image.
+
+procedure mg_setboundary2 (im, col1, col2, line1, line2, btype, bconst,
+	nxymargin)
+
+pointer	im			# IMIO pointer
+int	col1, col2		# Range of columns
+int	line1, line2		# Range of lines
+int	btype			# Boundary extension type
+real	bconst			# Constant for constant boundary extension
+int	nxymargin		# Number of margin pixels to allow
+
+int	btypes[5]
+int	nbndrypix
+
+data	btypes /BT_CONSTANT, BT_NEAREST, BT_REFLECT, BT_WRAP, BT_PROJECT/
+
+begin
+	nbndrypix = 0
+	nbndrypix = max (nbndrypix, 1 - col1)
+	nbndrypix = max (nbndrypix, col2 - IM_LEN(im, 1))
+	nbndrypix = max (nbndrypix, 1 - line1)
+	nbndrypix = max (nbndrypix, line2 - IM_LEN(im, 2))
+
+	call imseti (im, IM_TYBNDRY, btypes[btype])
+	call imseti (im, IM_NBNDRYPIX, nbndrypix + nxymargin + 1)
+	if (btypes[btype] == BT_CONSTANT)
+	    call imsetr (im, IM_BNDRYPIXVAL, bconst)
+end