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