aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/imgeom/src/t_magnify.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/imgeom/src/t_magnify.x')
-rw-r--r--pkg/images/imgeom/src/t_magnify.x624
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
generated by cgit (git 2.34.1) at 2025-09-20 05:18:10 -0400