aboutsummaryrefslogtreecommitdiff
path: root/noao/nproto/ir
diff options
context:
space:
mode:
authorJoseph Hunkeler <jhunkeler@gmail.com>2015-07-08 20:46:52 -0400
committerJoseph Hunkeler <jhunkeler@gmail.com>2015-07-08 20:46:52 -0400
commitfa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 (patch)
treebdda434976bc09c864f2e4fa6f16ba1952b1e555 /noao/nproto/ir
downloadiraf-linux-fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4.tar.gz
Initial commit
Diffstat (limited to 'noao/nproto/ir')
-rw-r--r--noao/nproto/ir/iralign.h55
-rw-r--r--noao/nproto/ir/iralign.x376
-rw-r--r--noao/nproto/ir/irdbio.x117
-rw-r--r--noao/nproto/ir/iriinit.x28
-rw-r--r--noao/nproto/ir/irimisec.x105
-rw-r--r--noao/nproto/ir/irimzero.x22
-rw-r--r--noao/nproto/ir/irindices.x139
-rw-r--r--noao/nproto/ir/irlinks.x496
-rw-r--r--noao/nproto/ir/irmatch1d.x122
-rw-r--r--noao/nproto/ir/irmatch2d.x276
-rw-r--r--noao/nproto/ir/irmedr.x35
-rw-r--r--noao/nproto/ir/iroverlap.x40
-rw-r--r--noao/nproto/ir/irqsort.x215
-rw-r--r--noao/nproto/ir/irtools.x147
-rw-r--r--noao/nproto/ir/mkpkg24
-rw-r--r--noao/nproto/ir/t_iralign.x134
-rw-r--r--noao/nproto/ir/t_irmatch1d.x159
-rw-r--r--noao/nproto/ir/t_irmatch2d.x159
-rw-r--r--noao/nproto/ir/t_irmosaic.x498
19 files changed, 3147 insertions, 0 deletions
diff --git a/noao/nproto/ir/iralign.h b/noao/nproto/ir/iralign.h
new file mode 100644
index 00000000..77cab3d4
--- /dev/null
+++ b/noao/nproto/ir/iralign.h
@@ -0,0 +1,55 @@
+# Header file for IR Mosaicing Routines
+
+# Define the structure
+
+define LEN_IRSTRUCT 35
+
+define IR_NCOLS Memi[$1] # x length of single subraster
+define IR_NROWS Memi[$1+1] # y length of a single subrasters
+define IR_NXOVERLAP Memi[$1+2] # x overlap between subrasters
+define IR_NYOVERLAP Memi[$1+3] # y overlap between subrasters
+define IR_NXSUB Memi[$1+4] # number of subrasters in x dimension
+define IR_NYSUB Memi[$1+5] # number of subrasters in y dimension
+define IR_NXRSUB Memi[$1+6] # x index of reference subraster
+define IR_NYRSUB Memi[$1+7] # y index of reference subraster
+define IR_XREF Memi[$1+8] # x offset of reference subraster
+define IR_YREF Memi[$1+9] # y offset of reference subraster
+define IR_CORNER Memi[$1+10] # starting corner for insertion
+define IR_ORDER Memi[$1+11] # row or column insertion
+define IR_RASTER Memi[$1+12] # raster order
+define IR_OVAL Memr[P2R($1+13)] # undefined value
+
+define IR_IC1 Memi[$1+14] # input image lower column limit
+define IR_IC2 Memi[$1+15] # input image upper column limit
+define IR_IL1 Memi[$1+16] # input image lower line limit
+define IR_IL2 Memi[$1+17] # input image upper line limit
+define IR_OC1 Memi[$1+18] # output image lower column limit
+define IR_OC2 Memi[$1+19] # output image upper column limit
+define IR_OL1 Memi[$1+20] # output image lower line limit
+define IR_OL2 Memi[$1+21] # output image upper line limit
+define IR_DELTAX Memi[$1+22] # x shifts
+define IR_DELTAY Memi[$1+23] # y shifts
+define IR_DELTAI Memi[$1+24] # intensity shifts
+
+define IR_XRSHIFTS Memi[$1+25] # x row links
+define IR_YRSHIFTS Memi[$1+26] # y row links
+define IR_NRSHIFTS Memi[$1+27] # number of row links
+define IR_XCSHIFTS Memi[$1+28] # x column links
+define IR_YCSHIFTS Memi[$1+29] # y column links
+define IR_NCSHIFTS Memi[$1+30] # number of column links
+
+# Define some useful constants
+
+define IR_LL 1
+define IR_LR 2
+define IR_UL 3
+define IR_UR 4
+
+define IR_ROW 1
+define IR_COLUMN 2
+
+define IR_COORDS 1
+define IR_SHIFTS 2
+define IR_FILE 3
+
+define MAX_NRANGES 100
diff --git a/noao/nproto/ir/iralign.x b/noao/nproto/ir/iralign.x
new file mode 100644
index 00000000..a1a431d5
--- /dev/null
+++ b/noao/nproto/ir/iralign.x
@@ -0,0 +1,376 @@
+include <imhdr.h>
+include "iralign.h"
+
+define NYOUT 16
+define NMARGIN 4
+
+# IR_SHIFTS -- Compute the input and output image column limits and the
+# x and y shifts.
+
+procedure ir_shifts (ir, im, outim, xrshifts, yrshifts, xcshifts,
+ ycshifts, ic1, ic2, il1, il2, oc1, oc2, ol1, ol2, deltax, deltay)
+
+pointer ir # pointer to the ir structure
+pointer im # pointer to the input image
+pointer outim # pointer to the output image
+real xrshifts[ARB] # x row shifts
+real yrshifts[ARB] # y row shifts
+real xcshifts[ARB] # x column shifts
+real ycshifts[ARB] # y column shifts
+int ic1[ARB] # input beginning column limits
+int ic2[ARB] # input ending column limits
+int il1[ARB] # input beginning line limits
+int il2[ARB] # input ending line limits
+int oc1[ARB] # output beginning column limits
+int oc2[ARB] # output ending column limits
+int ol1[ARB] # output beginning line limits
+int ol2[ARB] # output ending line limits
+real deltax[ARB] # x shifts
+real deltay[ARB] # x shifts
+
+
+int i, j, k, nimages, nxsize, nysize, nimcols, nimlines
+int c1ref, c2ref, l1ref, l2ref, ideltax, ideltay
+
+begin
+ # Find the position in the output image of the reference subraster.
+ nxsize = IR_NCOLS(ir) - IR_NXOVERLAP(ir)
+ nysize = IR_NROWS(ir) - IR_NYOVERLAP(ir)
+ c1ref = (IR_NXRSUB(ir) - 1) * nxsize + 1 + IR_XREF(ir)
+ c2ref = c1ref + IR_NCOLS(ir) - 1
+ l1ref = (IR_NYRSUB(ir) - 1) * nysize + 1 + IR_YREF(ir)
+ l2ref = l1ref + IR_NROWS(ir) - 1
+ nimages = IR_NXSUB(ir) * IR_NYSUB(ir)
+
+ # Extract the subrasters one by one.
+ do i = 1, nimages {
+
+ # Compute the indices of each subraster.
+ call ir_indices (i, j, k, IR_NXSUB(ir), IR_NYSUB(ir),
+ IR_CORNER(ir), IR_RASTER(ir), IR_ORDER(ir))
+
+ # Compute the indices of the input subraster.
+ nimcols = IM_LEN(im,1)
+ nimlines = IM_LEN(im,2)
+ ic1[i] = max (1, min (1 + (j - 1) * nxsize, nimcols))
+ ic2[i] = min (nimcols, max (1, ic1[i] + IR_NCOLS(ir) - 1))
+ il1[i] = max (1, min (1 + (k - 1) * nysize, nimlines))
+ il2[i] = min (nimlines, max (1, il1[i] + IR_NROWS(ir) - 1))
+
+ # Compute the shift relative to the input subraster.
+ call ir_mkshift (xrshifts, yrshifts, xcshifts, ycshifts,
+ IR_NXSUB(ir), IR_NYSUB(ir), j, k, IR_NXRSUB(ir),
+ IR_NYRSUB(ir), IR_ORDER(ir), deltax[i], deltay[i])
+ ideltax = nint (deltax[i])
+ ideltay = nint (deltay[i])
+
+ # Get the output buffer.
+ oc1[i] = c1ref + (j - IR_NXRSUB(ir)) * IR_NCOLS(ir) +
+ ideltax
+ oc2[i] = c2ref + (j - IR_NXRSUB(ir)) * IR_NCOLS(ir) +
+ ideltax
+ ol1[i] = l1ref + (k - IR_NYRSUB(ir)) * IR_NROWS(ir) +
+ ideltay
+ ol2[i] = l2ref + (k - IR_NYRSUB(ir)) * IR_NROWS(ir) +
+ ideltay
+ }
+end
+
+
+# IR_FSHIFTS -- Compute the input and output column limits.
+
+procedure ir_fshifts (ir, im, outim, deltax, deltay, ic1, ic2, il1, il2,
+ oc1, oc2, ol1, ol2)
+
+pointer ir # pointer to the ir structure
+pointer im # pointer to the input image
+pointer outim # pointer to the output image
+real deltax[ARB] # x shifts
+real deltay[ARB] # x shifts
+int ic1[ARB] # input beginning column limits
+int ic2[ARB] # input ending column limits
+int il1[ARB] # input beginning line limits
+int il2[ARB] # input ending line limits
+int oc1[ARB] # output beginning column limits
+int oc2[ARB] # output ending column limits
+int ol1[ARB] # output beginning line limits
+int ol2[ARB] # output ending line limits
+
+
+int i, j, k, nimages, nxsize, nysize, nimcols, nimlines
+int c1ref, c2ref, l1ref, l2ref, ideltax, ideltay
+
+begin
+ # Find the position in the output image of the reference subraster.
+ nxsize = IR_NCOLS(ir) - IR_NXOVERLAP(ir)
+ nysize = IR_NROWS(ir) - IR_NYOVERLAP(ir)
+ c1ref = (IR_NXRSUB(ir) - 1) * nxsize + 1 + IR_XREF(ir)
+ c2ref = c1ref + IR_NCOLS(ir) - 1
+ l1ref = (IR_NYRSUB(ir) - 1) * nysize + 1 + IR_YREF(ir)
+ l2ref = l1ref + IR_NROWS(ir) - 1
+ nimages = IR_NXSUB(ir) * IR_NYSUB(ir)
+
+ # Extract the subrasters one by one.
+ do i = 1, nimages {
+
+ # Compute the indices of each subraster.
+ call ir_indices (i, j, k, IR_NXSUB(ir), IR_NYSUB(ir),
+ IR_CORNER(ir), IR_RASTER(ir), IR_ORDER(ir))
+
+ # Compute the indices of the input subraster.
+ nimcols = IM_LEN(im,1)
+ nimlines = IM_LEN(im,2)
+ ic1[i] = max (1, min (1 + (j - 1) * nxsize, nimcols))
+ ic2[i] = min (nimcols, max (1, ic1[i] + IR_NCOLS(ir) - 1))
+ il1[i] = max (1, min (1 + (k - 1) * nysize, nimlines))
+ il2[i] = min (nimlines, max (1, il1[i] + IR_NROWS(ir) - 1))
+
+ # Compute the shift relative to the input subraster.
+ ideltax = nint (deltax[i])
+ ideltay = nint (deltay[i])
+
+ # Get the output buffer.
+ oc1[i] = c1ref + (j - IR_NXRSUB(ir)) * IR_NCOLS(ir) +
+ ideltax
+ oc2[i] = c2ref + (j - IR_NXRSUB(ir)) * IR_NCOLS(ir) +
+ ideltax
+ ol1[i] = l1ref + (k - IR_NYRSUB(ir)) * IR_NROWS(ir) +
+ ideltay
+ ol2[i] = l2ref + (k - IR_NYRSUB(ir)) * IR_NROWS(ir) +
+ ideltay
+ }
+end
+
+
+# IR_SUBALIGN -- Align all the subrasters.
+
+procedure ir_subalign (ir, im, outim, trimlimits, ic1, ic2, il1, il2,
+ oc1, oc2, ol1, ol2, deltax, deltay, deltai, match, interp, verbose)
+
+pointer ir # pointer to the ir structure
+pointer im # pointer to the input image
+pointer outim # pointer to the output image
+char trimlimits[ARB] # compute the trim section
+int ic1[ARB] # input image beginning columns
+int ic2[ARB] # input image ending columns
+int il1[ARB] # input image beginning rows
+int il2[ARB] # input image ending rows
+int oc1[ARB] # output image beginning columns
+int oc2[ARB] # output image ending columns
+int ol1[ARB] # output image beginning rows
+int ol2[ARB] # output image ending rows
+real deltax[ARB] # array of x shifts
+real deltay[ARB] # array of y shifts
+real deltai[ARB] # array of intensity shifts
+int match # match intensities ?
+int interp # type of interpolant
+int verbose # print messages
+
+int i, k, tl1, tl2, tc1, tc2, nimcols, nimlines, nimages
+int ideltax, ideltay, lxoffset, hxoffset, lyoffset, hyoffset
+int ixoffset, iyoffset, nocols, norows, cin1, cin2, nicols
+int tlin1, lin1, lin2, nilines, lout1, lout2, nyout, fstline, lstline
+pointer sp, x, y, msi, inbuf, outbuf, ptr
+real dx, dy, ytemp
+int ir_decode_section()
+pointer imps2r()
+
+begin
+ # Allocate temporary space.
+ call smark (sp)
+ call salloc (x, IR_NCOLS(ir), TY_REAL)
+ call salloc (y, IR_NCOLS(ir), TY_REAL)
+
+ # Decode the trimsection.
+ if (ir_decode_section (trimlimits, IR_NCOLS(ir), IR_NROWS(ir),
+ tc1, tc2, tl1, tl2) == ERR) {
+ tc1 = 0
+ tc2 = 0
+ tl1 = 0
+ tl2 = 0
+ } else {
+ tc1 = max (0, min (tc1, IR_NCOLS(ir)))
+ tc2 = max (0, min (tc2, IR_NCOLS(ir)))
+ tl1 = max (0, min (tl1, IR_NROWS(ir)))
+ tl2 = max (0, min (tl2, IR_NROWS(ir)))
+ }
+
+ # Initialize the interpolant.
+ call msiinit (msi, interp)
+
+ nimcols = IM_LEN(outim,1)
+ nimlines = IM_LEN(outim,2)
+
+ # Extract the subrasters one by one.
+ nimages = IR_NXSUB(ir) * IR_NYSUB(ir)
+ do i = 1, nimages {
+
+ inbuf = NULL
+
+ # Reject and subraster which is off the image.
+ if (oc1[i] > nimcols || oc2[i] < 1 || ol1[i] > nimlines ||
+ ol2[i] < 1)
+ next
+
+ # Compute the integer and fractional part of the shift.
+ ideltax = nint (deltax[i])
+ ideltay = nint (deltay[i])
+ dx = deltax[i] - ideltax
+ dy = deltay[i] - ideltay
+
+ # Compute the output image limits.
+ lxoffset = max (1 - oc1[i], tc1)
+ hxoffset = max (oc2[i] - nimcols, tc2)
+ oc1[i] = max (1, min (nimcols, oc1[i] + lxoffset))
+ oc2[i] = min (nimcols, max (1, oc2[i] - hxoffset))
+ nocols = oc2[i] - oc1[i] + 1
+ lyoffset = max (1 - ol1[i], tl1)
+ hyoffset = max (ol2[i] - nimlines, tl2)
+ ol1[i] = max (1, min (nimlines, ol1[i] + lyoffset))
+ ol2[i] = min (nimlines, max (1, ol2[i] - hyoffset))
+ norows = ol2[i] - ol1[i] + 1
+
+ # Compute some input image parameters.
+ cin1 = max (ic1[i], min (ic1[i] + lxoffset - NMARGIN, ic2[i]))
+ cin2 = min (ic2[i], max (ic2[i] - hxoffset + NMARGIN, ic1[i]))
+ nicols = cin2 - cin1 + 1
+
+ # Compute the x offset and x interpolation coordinates.
+ ixoffset = min (lxoffset, NMARGIN)
+ do k = 1, nicols
+ Memr[x+k-1] = max (1.0, min (real (nicols), real (k + ixoffset -
+ dx)))
+
+ # Subdivide the image and do the shifting.
+ for (lout1 = ol1[i]; lout1 <= ol2[i]; lout1 = lout1 + NYOUT) {
+
+ # Compute the output image limits.
+ lout2 = min (ol2[i], lout1 + NYOUT - 1)
+ nyout = lout2 - lout1 + 1
+
+ # Compute the input image limits.
+ tlin1 = il1[i] + lyoffset + lout1 - ol1[i]
+ lin2 = min (il2[i], max (tlin1 + nyout + NMARGIN - 1, il1[i]))
+ lin1 = max (il1[i], min (tlin1 - NMARGIN, il2[i]))
+ nilines = lin2 - lin1 + 1
+
+ # Get the appropriate input image section and fit the
+ # interpolant.
+ if ((inbuf == NULL) || (lin1 < fstline) || (lin2 > lstline)) {
+ fstline = lin1
+ lstline = lin2
+ call ir_buf (im, cin1, cin2, lin1, lin2, inbuf)
+ call msifit (msi, Memr[inbuf], nicols, nilines, nicols)
+ }
+
+ # Get the y offset and y interpolation coordinates.
+ #iyoffset = max (0, lout1 - ideltay - lin1)
+ if (lout1 == ol1[i])
+ iyoffset = min (lyoffset, NMARGIN)
+ else
+ iyoffset = tlin1 - lin1
+
+ # Shift the input images.
+ outbuf = imps2r (outim, oc1[i], oc2[i], lout1, lout2)
+ ptr = outbuf
+ do k = 1, nyout {
+ ytemp = max (1.0, min (real (nilines), real (k + iyoffset -
+ dy)))
+ call amovkr (ytemp, Memr[y], nocols)
+ call msivector (msi, Memr[x], Memr[y], Memr[ptr], nocols)
+ ptr = ptr + nocols
+ }
+
+ # Shift the intensities.
+ if (match == YES && ! IS_INDEFR(deltai[i]))
+ call aaddkr (Memr[outbuf], deltai[i], Memr[outbuf],
+ nocols * nyout)
+ }
+
+ if (inbuf != NULL)
+ call mfree (inbuf, TY_REAL)
+ inbuf = NULL
+
+ # Print a message.
+ if (verbose == YES) {
+ call printf (" %s[%d:%d,%d:%d] [%d:%d,%d:%d] %g %g")
+ call pargstr (IM_HDRFILE(im))
+ call pargi (ic1[i])
+ call pargi (ic2[i])
+ call pargi (il1[i])
+ call pargi (il2[i])
+ call pargi (lxoffset + 1)
+ call pargi (lxoffset + nocols)
+ call pargi (lyoffset + 1)
+ call pargi (lyoffset + norows)
+ call pargr (deltax[i])
+ call pargr (deltay[i])
+ call printf (" %s[%d:%d,%d:%d] %g\n")
+ call pargstr (IM_HDRFILE(outim))
+ call pargi (oc1[i])
+ call pargi (oc2[i])
+ call pargi (ol1[i])
+ call pargi (ol2[i])
+ call pargr (deltai[i])
+ }
+
+ }
+
+ call msifree (msi)
+ call sfree (sp)
+end
+
+
+# IR_BUF -- Procedure to provide a buffer of image lines with minimum reads.
+
+procedure ir_buf (im, col1, col2, line1, line2, buf)
+
+pointer im # pointer to input image
+int col1, col2 # column range of input buffer
+int line1, line2 # line range of input buffer
+pointer buf # buffer
+
+int i, ncols, nlines, nclast, llast1, llast2, nllast
+pointer buf1, buf2
+
+pointer imgs2r()
+
+begin
+ ncols = col2 - col1 + 1
+ nlines = line2 - line1 + 1
+
+ if (buf == NULL) {
+ call malloc (buf, ncols * nlines, TY_REAL)
+ llast1 = line1 - nlines
+ llast2 = line2 - nlines
+ } else if ((nlines != nllast) || (ncols != nclast)) {
+ call realloc (buf, ncols * nlines, TY_REAL)
+ llast1 = line1 - nlines
+ llast2 = line2 - nlines
+ }
+
+ if (line1 < llast1) {
+ do i = line2, line1, -1 {
+ if (i > llast1)
+ buf1 = buf + (i - llast1) * ncols
+ else
+ buf1 = imgs2r (im, col1, col2, i, i)
+ buf2 = buf + (i - line1) * ncols
+ call amovr (Memr[buf1], Memr[buf2], ncols)
+ }
+ } else if (line2 > llast2) {
+ do i = line1, line2 {
+ if (i < llast2)
+ buf1 = buf + (i - llast1) * ncols
+ else
+ buf1 = imgs2r (im, col1, col2, i, i)
+ buf2 = buf + (i - line1) * ncols
+ call amovr (Memr[buf1], Memr[buf2], ncols)
+ }
+ }
+
+ llast1 = line1
+ llast2 = line2
+ nclast = ncols
+ nllast = nlines
+end
diff --git a/noao/nproto/ir/irdbio.x b/noao/nproto/ir/irdbio.x
new file mode 100644
index 00000000..00e40532
--- /dev/null
+++ b/noao/nproto/ir/irdbio.x
@@ -0,0 +1,117 @@
+include "iralign.h"
+
+# IR_DTRPARAMS -- Procedure to read in the parameters from the database file.
+
+procedure ir_dtrparams (dt, image, ir)
+
+pointer dt # pointer to the database file
+char image[ARB] # input image
+pointer ir # pointer to the ir structure
+
+int recnum, nsubrasters
+pointer sp, str
+int dtlocate(), dtgeti(), strmatch()
+real dtgetr()
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+ recnum = dtlocate (dt, image)
+
+ IR_NCOLS(ir) = dtgeti (dt, recnum, "ncols")
+ IR_NROWS(ir) = dtgeti (dt, recnum, "nrows")
+ IR_NXSUB(ir) = dtgeti (dt, recnum, "nxsub")
+ IR_NYSUB(ir) = dtgeti (dt, recnum, "nysub")
+ IR_NXOVERLAP(ir) = dtgeti (dt, recnum, "nxoverlap")
+ IR_NYOVERLAP(ir) = dtgeti (dt, recnum, "nyoverlap")
+
+ call dtgstr (dt, recnum, "corner", Memc[str], SZ_FNAME)
+ if (strmatch (Memc[str], "ll") != 0)
+ IR_CORNER(ir) = IR_LL
+ else if (strmatch (Memc[str], "lr") != 0)
+ IR_CORNER(ir) = IR_LR
+ else if (strmatch (Memc[str], "ul") != 0)
+ IR_CORNER(ir) = IR_UL
+ else if (strmatch (Memc[str], "ur") != 0)
+ IR_CORNER(ir) = IR_UR
+ else
+ IR_CORNER(ir) = IR_LL
+
+ call dtgstr (dt, recnum, "order", Memc[str], SZ_FNAME)
+ if (strmatch (Memc[str], "column") != 0)
+ IR_ORDER(ir) = IR_COLUMN
+ else if (strmatch (Memc[str], "row") != 0)
+ IR_ORDER(ir) = IR_ROW
+ else
+ IR_ORDER(ir) = IR_ROW
+
+ call dtgstr (dt, recnum, "raster", Memc[str], SZ_FNAME)
+ if (strmatch (Memc[str], "yes") != 0)
+ IR_RASTER(ir) = YES
+ else if (strmatch (Memc[str], "no") != 0)
+ IR_RASTER(ir) = NO
+ else
+ IR_RASTER(ir) = NO
+
+ IR_OVAL(ir) = dtgetr (dt, recnum, "oval")
+ nsubrasters = dtgeti (dt, recnum, "nsubrasters")
+
+ call sfree (sp)
+end
+
+
+# IR_DTWPARAMS -- Procedure to write out the parameters to the output file
+
+procedure ir_dtwparams (dt, outimage, trimsection, medsection, ir)
+
+pointer dt # pointer to the database file
+char outimage[ARB] # name of the output image
+char trimsection[ARB]# input subraster section
+char medsection[ARB] # section for computing the median
+pointer ir # pointer to the ir structure
+
+bool itob()
+
+begin
+ call dtptime (dt)
+ call dtput (dt, "begin\t%s\n")
+ call pargstr (outimage)
+ call dtput (dt, "\ttrimsection\t%s\n")
+ call pargstr (trimsection)
+ call dtput (dt, "\tmedsection\t\t%s\n")
+ call pargstr (medsection)
+ call dtput (dt, "\tncols\t\t%d\n")
+ call pargi (IR_NCOLS(ir))
+ call dtput (dt, "\tnrows\t\t%d\n")
+ call pargi (IR_NROWS(ir))
+ call dtput (dt, "\tnxsub\t\t%d\n")
+ call pargi (IR_NXSUB(ir))
+ call dtput (dt, "\tnysub\t\t%d\n")
+ call pargi (IR_NYSUB(ir))
+ call dtput (dt, "\tnxoverlap\t%d\n")
+ call pargi (IR_NXOVERLAP(ir))
+ call dtput (dt, "\tnyoverlap\t%d\n")
+ call pargi (IR_NYOVERLAP(ir))
+ call dtput (dt, "\tcorner\t\t%s\n")
+ switch (IR_CORNER(ir)) {
+ case IR_LL:
+ call pargstr ("ll")
+ case IR_LR:
+ call pargstr ("lr")
+ case IR_UL:
+ call pargstr ("ul")
+ case IR_UR:
+ call pargstr ("ur")
+ }
+ call dtput (dt, "\torder\t\t%s\n")
+ switch (IR_ORDER(ir)) {
+ case IR_ROW:
+ call pargstr ("row")
+ case IR_COLUMN:
+ call pargstr ("column")
+ }
+ call dtput (dt, "\traster\t\t%b\n")
+ call pargb (itob (IR_RASTER(ir)))
+ call dtput (dt, "\toval\t\t%g\n")
+ call pargr (IR_OVAL(ir))
+end
diff --git a/noao/nproto/ir/iriinit.x b/noao/nproto/ir/iriinit.x
new file mode 100644
index 00000000..a97ade8e
--- /dev/null
+++ b/noao/nproto/ir/iriinit.x
@@ -0,0 +1,28 @@
+# IR_VECINIT -- Procedure to initialize the intensity matching algorithm.
+# If the ranges are undefined and no matching is to take place the
+# ishifts are set to INDEFR and the routine returns. Otherwise the shifts
+# are all initialized to zero and shifts for the missing subrasters are
+# set to INDEFR.
+
+procedure ir_vecinit (deltai, nsubrasters, ranges)
+
+real deltai[ARB] # intensity shifts
+int nsubrasters # number of subrasters
+int ranges[ARB] # ranges of missing subrasters
+
+int num
+int get_next_number()
+
+begin
+ # Initialize the shifts to INDEFR.
+ call amovkr (INDEFR, deltai, nsubrasters)
+ if (ranges[1] == NULL)
+ return
+
+ num = 0
+ while (get_next_number (ranges, num) != EOF) {
+ if (num > nsubrasters)
+ break
+ deltai[num] = 0.0
+ }
+end
diff --git a/noao/nproto/ir/irimisec.x b/noao/nproto/ir/irimisec.x
new file mode 100644
index 00000000..1b6936db
--- /dev/null
+++ b/noao/nproto/ir/irimisec.x
@@ -0,0 +1,105 @@
+include <ctype.h>
+
+# IR_DECODE_SECTION -- Procedure to decode the reference section.
+
+int procedure ir_decode_section (section, ncols, nrows, c1ref, c2ref, l1ref,
+ l2ref)
+
+char section[ARB] # reference subraster section
+int ncols # number of columns in the image
+int nrows # number of rows in the image
+int c1ref # initial column
+int c2ref # final reference column
+int l1ref # initial reference line
+int l2ref # final reference line
+
+char leftbkt
+int index, ip, step
+int ir_decode_subscript(), stridx()
+
+begin
+ leftbkt = '['
+ index = stridx (leftbkt, section)
+ if (index == 0)
+ return (ERR)
+ ip = index + 1
+ if (ir_decode_subscript (section, ip, ncols, c1ref, c2ref, step) == ERR)
+ return (ERR)
+ if (ir_decode_subscript (section, ip, nrows, l1ref, l2ref, step) == ERR)
+ return (ERR)
+ return (OK)
+end
+
+
+# IR_DECODE_SUBSCRIPT -- Decode a single subscript expression to produce the
+# range of values for that subscript (X1:X2), and the sampling step size, STEP.
+# Note that X1 may be less than, greater than, or equal to X2, and STEP may
+# be a positive or negative nonzero integer. Various shorthand notations are
+# permitted, as is embedded whitespace.
+
+int procedure ir_decode_subscript (section, ip, maxnumber, x1, x2, step)
+
+char section[ARB]
+int ip
+int maxnumber
+long x1, x2, step, temp
+int ctol()
+
+begin
+ x1 = 1
+ x2 = maxnumber
+ step = 1
+
+ while (IS_WHITE(section[ip]))
+ ip = ip + 1
+
+ # Get X1, X2.
+ if (ctol (section, ip, temp) > 0) { # [x1
+ x1 = temp
+ if (section[ip] == ':') {
+ ip = ip + 1
+ if (ctol (section, ip, x2) == 0) # [x1:x2
+ return (ERR)
+ } else
+ x2 = x1
+
+ } else if (section[ip] == '-') {
+ x1 = maxnumber # [-*
+ x2 = 1
+ ip = ip + 1
+ if (section[ip] == '*')
+ ip = ip + 1
+
+ } else if (section[ip] == '*') # [*
+ ip = ip + 1
+
+ while (IS_WHITE(section[ip]))
+ ip = ip + 1
+
+ # Get sample step size, if give.
+ if (section[ip] == ':') { # ..:step
+ ip = ip + 1
+ if (ctol (section, ip, step) == 0)
+ return (ERR)
+ else if (step == 0)
+ return (ERR)
+ }
+
+ # Allow notation such as "-*:5", (or even "-:5") where the step
+ # is obviously supposed to be negative.
+
+ if (x1 > x2 && step > 0)
+ step = -step
+
+ while (IS_WHITE(section[ip]))
+ ip = ip + 1
+
+ if (section[ip] == ',') {
+ ip = ip + 1
+ return (OK)
+ } else if (section[ip] == ']')
+ return (OK)
+ else
+ return (ERR)
+
+end
diff --git a/noao/nproto/ir/irimzero.x b/noao/nproto/ir/irimzero.x
new file mode 100644
index 00000000..013ab2d6
--- /dev/null
+++ b/noao/nproto/ir/irimzero.x
@@ -0,0 +1,22 @@
+
+# IR_IMZERO -- Fill the output image with a constant value.
+
+procedure ir_imzero (im, ncols, nlines, value)
+
+pointer im # pointer to the output image
+int ncols # number of columns
+int nlines # number of lines
+real value # default blank value
+
+int i
+pointer obuf
+pointer impl2r()
+
+begin
+ do i = 1, nlines {
+ obuf = impl2r (im, i)
+ if (obuf == EOF)
+ call error (0, "Error writing output image.")
+ call amovkr (value, Memr[obuf], ncols)
+ }
+end
diff --git a/noao/nproto/ir/irindices.x b/noao/nproto/ir/irindices.x
new file mode 100644
index 00000000..2ff94d81
--- /dev/null
+++ b/noao/nproto/ir/irindices.x
@@ -0,0 +1,139 @@
+include "iralign.h"
+
+# IR_INDICES -- Given the number in the list for a missing subraster and
+# information about how the subrasters were written return the i and j
+# indices of the specified subrasters.
+
+procedure ir_indices (num, i, j, nxsub, nysub, corner, raster, order)
+
+int num # number of the subraster
+int i,j # indices of the subraster
+int nxsub,nysub # number of subrasters in x and y
+int corner # starting corner
+int raster # raster order
+int order # column or row order
+
+begin
+ switch (corner) {
+ case IR_LL:
+ if (order == IR_ROW) {
+ if (mod (num, nxsub) == 0) {
+ j = num / nxsub
+ if (raster == YES && mod (j,2) == 0)
+ i = 1
+ else
+ i = nxsub
+ } else {
+ j = num / nxsub + 1
+ if (raster == YES && mod (j,2) == 0)
+ i = nxsub - mod (num, nxsub) + 1
+ else
+ i = mod (num, nxsub)
+ }
+ } else if (order == IR_COLUMN) {
+ if (mod (num, nysub) == 0) {
+ i = num / nysub
+ if (raster == YES && mod (i,2) == 0)
+ j = 1
+ else
+ j = nysub
+ } else {
+ i = num / nysub + 1
+ if (raster == YES && mod (i,2) == 0)
+ j = nysub - mod (num, nysub) + 1
+ else
+ j = mod (num, nysub)
+ }
+ }
+ case IR_LR:
+ if (order == IR_ROW) {
+ if (mod (num, nxsub) == 0) {
+ j = num / nxsub
+ if (raster == YES && mod (j,2) == 0)
+ i = nxsub
+ else
+ i = 1
+ } else {
+ j = num / nxsub + 1
+ if (raster == YES && mod (j,2) == 0)
+ i = mod (num, nxsub)
+ else
+ i = nxsub - mod (num, nxsub) + 1
+ }
+ } else if (order == IR_COLUMN) {
+ if (mod (num, nysub) == 0) {
+ i = nxsub - num / nysub + 1
+ if (raster == YES && mod (i,2) != 0)
+ j = 1
+ else
+ j = nysub
+ } else {
+ i = nxsub - num / nysub
+ if (raster == YES && mod (i,2) != 0)
+ j = nysub - mod (num, nysub) + 1
+ else
+ j = mod (num, nysub)
+ }
+ }
+ case IR_UL:
+ if (order == IR_ROW) {
+ if (mod (num, nxsub) == 0) {
+ j = nysub - num / nxsub + 1
+ if (raster == YES && mod (j,2) != 0)
+ i = 1
+ else
+ i = nxsub
+ } else {
+ j = nysub - num / nxsub
+ if (raster == YES && mod (j,2) != 0)
+ i = nxsub - mod (num, nxsub) + 1
+ else
+ i = mod (num, nxsub)
+ }
+ } else if (order == IR_COLUMN) {
+ if (mod (num, nysub) == 0) {
+ i = num / nysub
+ if (raster == YES && mod (i,2) == 0)
+ j = nysub
+ else
+ j = 1
+ } else {
+ i = num / nysub + 1
+ if (raster == YES && mod (i,2) == 0)
+ j = mod (num, nysub)
+ else
+ j = nysub - mod (num, nysub) + 1
+ }
+ }
+ case IR_UR:
+ if (order == IR_ROW) {
+ if (mod (num, nxsub) == 0) {
+ j = nysub - num / nxsub + 1
+ if (raster == YES && mod (j,2) != 0)
+ i = nxsub
+ else
+ i = 1
+ } else {
+ j = nysub - num / nxsub
+ if (raster == YES && mod (j,2) != 0)
+ i = mod (num, nxsub)
+ else
+ i = nxsub - mod (num, nxsub) + 1
+ }
+ } else if (order == IR_COLUMN) {
+ if (mod (num, nysub) == 0) {
+ i = nxsub - num / nysub + 1
+ if (raster == YES && mod (i,2) != 0)
+ j = nysub
+ else
+ j = 1
+ } else {
+ i = nxsub - num / nysub
+ if (raster == YES && mod (i,2) != 0)
+ j = mod (num, nysub)
+ else
+ j = nysub - mod (num, nysub) + 1
+ }
+ }
+ }
+end
diff --git a/noao/nproto/ir/irlinks.x b/noao/nproto/ir/irlinks.x
new file mode 100644
index 00000000..2b8b3a4a
--- /dev/null
+++ b/noao/nproto/ir/irlinks.x
@@ -0,0 +1,496 @@
+include "iralign.h"
+
+# IR_LINKS -- Procedure to compute the shifts for each subraster.
+
+int procedure ir_links (cl, xrshift, yrshift, xcshift, ycshift, nrshift,
+ ncshift, ncols, nrows, nxrsub, nyrsub, nxsub, nysub, nxoverlap,
+ nyoverlap, order)
+
+int cl # coordinate list descriptor
+real xrshift[nxsub,ARB] # x row shifts
+real yrshift[nxsub,ARB] # y row shifts
+real xcshift[nxsub,ARB] # x column shifts
+real ycshift[nxsub,ARB] # y column shifts
+int nrshift[nxsub,ARB] # number of row shifts
+int ncshift[nxsub,ARB] # number of column shifts
+int ncols # number of columns per subraster
+int nrows # number of rows per subraster
+int nxrsub # column index of reference subraster
+int nyrsub # row index of reference subraster
+int nxsub # number of subrasters in x
+int nysub # number of subrasters in y
+int nxoverlap # number of columns of overlap
+int nyoverlap # number of rows of overlap
+int order # row or column order
+
+int i, j, nxsize, nysize, ilimit, olimit, nshifts
+pointer sp, xcolavg, ycolavg, xrowavg, yrowavg, nrowavg, ncolavg
+real isign, jsign, xrmed, yrmed, xcmed, ycmed
+int ir_decode_shifts()
+real irmedr()
+
+begin
+ # Allocate temporary space.
+ if (order == IR_COLUMN) {
+ ilimit = nysub
+ olimit = nxsub
+ } else {
+ ilimit = nxsub
+ olimit = nysub
+ }
+
+ # Clear the shift arrays.
+ call aclrr (xrshift, nxsub * nysub)
+ call aclrr (yrshift, nxsub * nysub)
+ call aclrr (xcshift, nxsub * nysub)
+ call aclrr (ycshift, nxsub * nysub)
+ call aclri (nrshift, nxsub * nysub)
+ call aclri (ncshift, nxsub * nysub)
+
+ # Accumulate the shifts.
+ nxsize = ncols - nxoverlap
+ nysize = nrows - nyoverlap
+ nshifts = ir_decode_shifts (cl, xrshift, yrshift, nrshift, xcshift,
+ ycshift, ncshift, nxsub, nysub, nxrsub, nyrsub, nxoverlap,
+ nyoverlap, nxsize, nysize)
+ if (nshifts == 0)
+ return (0)
+
+ # Allocate working space.
+ call smark (sp)
+ call salloc (xcolavg, olimit, TY_REAL)
+ call salloc (ycolavg, olimit, TY_REAL)
+ call salloc (ncolavg, olimit, TY_INT)
+ call salloc (xrowavg, olimit, TY_REAL)
+ call salloc (yrowavg, olimit, TY_REAL)
+ call salloc (nrowavg, olimit, TY_INT)
+
+ # Clear the accumulator arrays.
+ call aclrr (Memr[xcolavg], olimit)
+ call aclrr (Memr[ycolavg], olimit)
+ call aclri (Memi[ncolavg], olimit)
+ call aclrr (Memr[xrowavg], olimit)
+ call aclrr (Memr[yrowavg], olimit)
+ call aclri (Memi[nrowavg], olimit)
+
+ # Compute the row or column sums.
+ if (order == IR_COLUMN) {
+ do i = 1, nxsub {
+ do j = 1, nysub {
+ if (nrshift[i,j] > 0) {
+ Memr[xrowavg+i-1] = Memr[xrowavg+i-1] +
+ abs (xrshift[i,j])
+ Memr[yrowavg+i-1] = Memr[yrowavg+i-1] +
+ abs (yrshift[i,j])
+ Memi[nrowavg+i-1] = Memi[nrowavg+i-1] + 1
+ }
+ if (ncshift[i,j] > 0) {
+ Memr[xcolavg+i-1] = Memr[xcolavg+i-1] +
+ abs (xcshift[i,j])
+ Memr[ycolavg+i-1] = Memr[ycolavg+i-1] +
+ abs (ycshift[i,j])
+ Memi[ncolavg+i-1] = Memi[ncolavg+i-1] + 1
+ }
+ }
+ }
+ } else {
+ do i = 1, nysub {
+ do j = 1, nxsub {
+ if (nrshift[j,i] > 0) {
+ Memr[xrowavg+i-1] = Memr[xrowavg+i-1] +
+ abs (xrshift[j,i])
+ Memr[yrowavg+i-1] = Memr[yrowavg+i-1] +
+ abs (yrshift[j,i])
+ Memi[nrowavg+i-1] = Memi[nrowavg+i-1] + 1
+ }
+ if (ncshift[j,i] > 0) {
+ Memr[xcolavg+i-1] = Memr[xcolavg+i-1] +
+ abs (xcshift[j,i])
+ Memr[ycolavg+i-1] = Memr[ycolavg+i-1] +
+ abs (ycshift[j,i])
+ Memi[ncolavg+i-1] = Memi[ncolavg+i-1] + 1
+ }
+ }
+ }
+ }
+
+ # Compute the averages.
+ do i = 1, olimit {
+ if (Memi[nrowavg+i-1] > 0) {
+ Memr[xrowavg+i-1] = Memr[xrowavg+i-1] / Memi[nrowavg+i-1]
+ Memr[yrowavg+i-1] = Memr[yrowavg+i-1] / Memi[nrowavg+i-1]
+ }
+ if (Memi[ncolavg+i-1] > 0) {
+ Memr[xcolavg+i-1] = Memr[xcolavg+i-1] / Memi[ncolavg+i-1]
+ Memr[ycolavg+i-1] = Memr[ycolavg+i-1] / Memi[ncolavg+i-1]
+ }
+ }
+
+ # Compute the medians of the row and column averages.
+ xrmed = irmedr (Memr[xrowavg], Memi[nrowavg], olimit)
+ yrmed = irmedr (Memr[yrowavg], Memi[nrowavg], olimit)
+ xcmed = irmedr (Memr[xcolavg], Memi[ncolavg], olimit)
+ ycmed = irmedr (Memr[ycolavg], Memi[ncolavg], olimit)
+
+ # Use the average shifts for subrasters with no information.
+ do j = 1, nysub {
+
+ if (j == nyrsub)
+ jsign = 0.0
+ else if (j < nyrsub)
+ jsign = 1.0
+ else
+ jsign = -1.0
+
+ do i = 1, nxsub {
+
+ if (i == nxrsub)
+ isign = 0.0
+ else if (i < nxrsub)
+ isign = 1.0
+ else
+ isign = -1.0
+
+ if (nrshift[i,j] <= 0) {
+ if (Memi[nrowavg+i-1] <= 0) {
+ xrshift[i,j] = isign * xrmed
+ yrshift[i,j] = jsign * yrmed
+ } else if (order == IR_COLUMN) {
+ xrshift[i,j] = isign * Memr[xrowavg+i-1]
+ yrshift[i,j] = jsign * Memr[yrowavg+i-1]
+ } else {
+ xrshift[i,j] = isign * Memr[xrowavg+j-1]
+ yrshift[i,j] = jsign * Memr[yrowavg+j-1]
+ }
+ }
+
+ if (ncshift[i,j] <= 0) {
+ if (Memi[ncolavg+i-1] <= 0) {
+ xcshift[i,j] = isign * xcmed
+ ycshift[i,j] = jsign * ycmed
+ } else if (order == IR_COLUMN) {
+ xcshift[i,j] = isign * Memr[xcolavg+i-1]
+ ycshift[i,j] = jsign * Memr[ycolavg+i-1]
+ } else {
+ xcshift[i,j] = isign * Memr[xcolavg+j-1]
+ ycshift[i,j] = jsign * Memr[ycolavg+j-1]
+ }
+ }
+ }
+ }
+
+ call sfree (sp)
+ return (nshifts)
+end
+
+
+# IR_DECODE_SHIFTS -- Procedure to accumulate shifts for each subraster.
+
+int procedure ir_decode_shifts (cl, xrshift, yrshift, nrshift, xcshift,
+ ycshift, ncshift, nxsub, nysub, nxrsub, nyrsub, nxoverlap,
+ nyoverlap, nxsize, nysize)
+
+int cl # coordinate list descriptor
+real xrshift[nxsub,ARB] # x row shifts
+real yrshift[nxsub,ARB] # y row shifts
+int nrshift[nxsub,ARB] # number of row shifts
+real xcshift[nxsub,ARB] # x column shifts
+real ycshift[nxsub,ARB] # y column shifts
+int ncshift[nxsub,ARB] # number of column shifts
+int nxsub # number of subrasters in x
+int nysub # number of subrasters in y
+int nxrsub # column index of reference subraster
+int nyrsub # row index of reference subraster
+int nxoverlap # number of columns of overlap
+int nyoverlap # number of rows of overlap
+int nxsize # size of unoverlapped region
+int nysize # size of unoverlapped region
+
+int i, j, nx1, ny1, nx2, ny2, r21, r22, stat, nshifts
+real x1, y1, x2, y2, xdif, xdifm, ydif, ydifm
+int fscan(), nscan()
+
+begin
+ nshifts = 0
+ while (fscan (cl) != EOF) {
+
+ # Get the first coordinate pair.
+ call gargr (x1)
+ call gargr (y1)
+ if (nscan () != 2)
+ next
+
+ # Compute which subraster 1 belongs to.
+ if (mod (int (x1), nxsize) == 0)
+ nx1 = int (x1) / nxsize
+ else
+ nx1 = int (x1) / nxsize + 1
+
+ if (mod (int (y1), nysize) == 0)
+ ny1 = int (y1) / nysize
+ else
+ ny1 = int (y1) / nysize + 1
+
+ # Get the second coordinate pair.
+ repeat {
+
+ stat = fscan (cl)
+ if (stat == EOF)
+ break
+ call gargr (x2)
+ call gargr (y2)
+
+ # Compute which subraster 2 belongs to.
+ if (nscan () == 2) {
+ if (mod (int (x2), nxsize) == 0)
+ nx2 = int (x2) / nxsize
+ else
+ nx2 = int (x2) / nxsize + 1
+ if (mod (int (y2), nysize) == 0)
+ ny2 = int (y2) / nysize
+ else
+ ny2 = int (y2) / nysize + 1
+ }
+
+ } until (nscan () == 2)
+ if (stat == EOF || nscan() != 2)
+ break
+
+ r21 = (nx1 - nxrsub) ** 2 + (ny1 - nyrsub) ** 2
+ r22 = (nx2 - nxrsub) ** 2 + (ny2 - nyrsub) ** 2
+
+ # Illegal shift
+ if (r21 == r22)
+ next
+
+ # Compute the shift for the first subraster.
+ else if (r21 > r22) {
+
+ xdif = x2 - x1
+ if (nxoverlap < 0) {
+ if (xdif < 0.0)
+ xdifm = xdif - nxoverlap
+ else if (xdif > 0.0)
+ xdifm = xdif + nxoverlap
+ } else
+ xdifm = xdif
+
+ ydif = y2 - y1
+ if (nyoverlap < 0) {
+ if (ydif < 0.0)
+ ydifm = ydif - nyoverlap
+ else if (ydif > 0.0)
+ ydifm = ydif + nyoverlap
+ } else
+ ydifm = ydif
+
+ if (nx1 == nx2) {
+ xcshift[nx1,ny1] = xcshift[nx1,ny1] + xdif
+ ycshift[nx1,ny1] = ycshift[nx1,ny1] + ydifm
+ ncshift[nx1,ny1] = ncshift[nx1,ny1] + 1
+ } else if (ny1 == ny2) {
+ xrshift[nx1,ny1] = xrshift[nx1,ny1] + xdifm
+ yrshift[nx1,ny1] = yrshift[nx1,ny1] + ydif
+ nrshift[nx1,ny1] = nrshift[nx1,ny1] + 1
+ } else
+ next
+
+ # Compute the shift for the second subraster.
+ } else {
+
+ xdif = x1 - x2
+ if (nxoverlap < 0) {
+ if (xdif < 0.0)
+ xdifm = xdif - nxoverlap
+ else if (xdif > 0.0)
+ xdifm = xdif + nxoverlap
+ } else
+ xdifm = xdif
+
+ ydif = y1 - y2
+ if (nyoverlap < 0) {
+ if (ydif < 0.0)
+ ydifm = ydif - nyoverlap
+ else if (ydif > 0.0)
+ ydifm = ydif + nyoverlap
+ } else
+ ydifm = ydif
+
+ if (nx1 == nx2) {
+ xcshift[nx2,ny2] = xcshift[nx2,ny2] + xdif
+ ycshift[nx2,ny2] = ycshift[nx2,ny2] + ydifm
+ ncshift[nx2,ny2] = ncshift[nx2,ny2] + 1
+ } else if (ny1 == ny2) {
+ xrshift[nx2,ny2] = xrshift[nx2,ny2] + xdifm
+ yrshift[nx2,ny2] = yrshift[nx2,ny2] + ydif
+ nrshift[nx2,ny2] = nrshift[nx2,ny2] + 1
+ } else
+ next
+ }
+
+ nshifts = nshifts + 1
+ }
+
+ # Compute the final shifts.
+ do j = 1, nysub {
+ do i = 1, nxsub {
+ if (nrshift[i,j] > 0) {
+ xrshift[i,j] = xrshift[i,j] / nrshift[i,j]
+ yrshift[i,j] = yrshift[i,j] / nrshift[i,j]
+ }
+ if (ncshift[i,j] > 0) {
+ xcshift[i,j] = xcshift[i,j] / ncshift[i,j]
+ ycshift[i,j] = ycshift[i,j] / ncshift[i,j]
+ }
+ }
+ }
+
+ return (nshifts)
+end
+
+
+# IR_CLINKS -- Procedure to compute the shifts for each subraster.
+
+int procedure ir_clinks (xrshift, yrshift, xcshift, ycshift, nxrsub, nyrsub,
+ nxsub, nysub, xshift, yshift)
+
+real xrshift[nxsub,ARB] # x row shifts
+real yrshift[nxsub,ARB] # y row shifts
+real xcshift[nxsub,ARB] # x column shifts
+real ycshift[nxsub,ARB] # y column shifts
+int nxrsub # x index of reference subraster
+int nyrsub # y index of reference subraster
+int nxsub # number of subrasters in x direction
+int nysub # number of subrasters in y direction
+real xshift # xshift of the coordinates
+real yshift # yshift of the coordinates
+
+int i, j, isign, jsign
+
+begin
+ do j = 1, nysub {
+ if (j == nyrsub)
+ jsign = 0
+ else if (j < nyrsub)
+ jsign = 1
+ else
+ jsign = -1
+
+ do i = 1, nxsub {
+ if (i == nxrsub)
+ isign = 0
+ else if (i < nxrsub)
+ isign = 1
+ else
+ isign = -1
+
+ xrshift[i,j] = isign * abs (xshift)
+ yrshift[i,j] = 0.0
+ xcshift[i,j] = 0.0
+ ycshift[i,j] = jsign * abs (yshift)
+ }
+ }
+
+ return (1)
+end
+
+
+# IR_FLINKS -- Routine to fetch the shifts directly
+
+int procedure ir_flinks (cl, deltax, deltay, deltai, max_nshifts)
+
+int cl # shifts file descriptor
+real deltax[ARB] # x shifts
+real deltay[ARB] # y shifts
+real deltai[ARB] # intensity shifts
+int max_nshifts # maximum number of shifts
+
+int nshifts
+int fscan(), nscan()
+
+begin
+ nshifts = 0
+ while ((fscan (cl) != EOF) && (nshifts < max_nshifts)) {
+ call gargr (deltax[nshifts+1])
+ call gargr (deltay[nshifts+1])
+ call gargr (deltai[nshifts+1])
+ if (nscan() < 2)
+ next
+ if (nscan() < 3)
+ deltai[nshifts+1] = 0.0
+ nshifts = nshifts + 1
+ }
+
+ return (nshifts)
+end
+
+
+# IR_MKSHIFT -- Routine to compute the total shift for each subraster.
+
+procedure ir_mkshift (xrshift, yrshift, xcshift, ycshift, nxsub, nysub,
+ xsubindex, ysubindex, nxrsub, nyrsub, order, deltax, deltay)
+
+real xrshift[nxsub,ARB] # x row shifts
+real yrshift[nxsub,ARB] # y row shifts
+real xcshift[nxsub,ARB] # x column shifts
+real ycshift[nxsub,ARB] # y column shifts
+int nxsub # number of subrasters in x direction
+int nysub # number of subrasters in y direction
+int xsubindex # x index of the subraster
+int ysubindex # y index of the subraster
+int nxrsub # x index of reference subraster
+int nyrsub # y index of reference subraster
+int order # row or column order
+real deltax # total x shift
+real deltay # total y shift
+
+int j
+
+begin
+ deltax = 0.0
+ deltay = 0.0
+
+ if (order == IR_COLUMN) {
+ if (ysubindex < nyrsub)
+ do j = ysubindex, nyrsub - 1 {
+ deltax = deltax + xcshift[xsubindex,j]
+ deltay = deltay + ycshift[xsubindex,j]
+ }
+ else if (ysubindex > nyrsub)
+ do j = nyrsub + 1, ysubindex {
+ deltax = deltax + xcshift[xsubindex,j]
+ deltay = deltay + ycshift[xsubindex,j]
+ }
+ if (xsubindex < nxrsub)
+ do j = xsubindex, nxrsub - 1 {
+ deltax = deltax + xrshift[j,nyrsub]
+ deltay = deltay + yrshift[j,nyrsub]
+ }
+ else if (xsubindex > nxrsub)
+ do j = nxrsub + 1, xsubindex {
+ deltax = deltax + xrshift[j,nyrsub]
+ deltay = deltay + yrshift[j,nyrsub]
+ }
+ } else {
+ if (xsubindex < nxrsub)
+ do j = xsubindex, nxrsub - 1{
+ deltax = deltax + xrshift[j,ysubindex]
+ deltay = deltay + yrshift[j,ysubindex]
+ }
+ else if (xsubindex > nxrsub)
+ do j = nxrsub + 1, xsubindex {
+ deltax = deltax + xrshift[j,ysubindex]
+ deltay = deltay + yrshift[j,ysubindex]
+ }
+ if (ysubindex < nyrsub)
+ do j = ysubindex, nyrsub - 1 {
+ deltax = deltax + xcshift[nxrsub,j]
+ deltay = deltay + ycshift[nxrsub,j]
+ }
+ else if (ysubindex > nyrsub)
+ do j = nyrsub + 1, ysubindex {
+ deltax = deltax + xcshift[nxrsub,j]
+ deltay = deltay + ycshift[nxrsub,j]
+ }
+ }
+end
diff --git a/noao/nproto/ir/irmatch1d.x b/noao/nproto/ir/irmatch1d.x
new file mode 100644
index 00000000..b3c6cdfb
--- /dev/null
+++ b/noao/nproto/ir/irmatch1d.x
@@ -0,0 +1,122 @@
+include <imhdr.h>
+include <pkg/dttext.h>
+include "iralign.h"
+
+# IR_M1MATCH -- Procedure to match images in the direction of observation
+# direction.
+
+procedure ir_m1match (ir, im, ranges, ic1, ic2, il1, il2, deltax, deltay,
+ deltai)
+
+pointer ir # pointer to the ir strucuture
+pointer im # pointer to the input image
+int ranges[ARB] # array elements to be skipped
+int ic1[ARB] # input beginning column limits
+int ic2[ARB] # output beginning column limits
+int il1[ARB] # input beginning line limits
+int il2[ARB] # output beginning line limits
+real deltax[ARB] # x shifts
+real deltay[ARB] # y shifts
+real deltai[ARB] # intensity shifts
+
+int num, nmod, turn_corner
+int pc1, pc2, pl1, pl2, c1, c2, l1, l2
+int pideltax, pideltay, ideltax, ideltay
+int oc1, oc2, ol1, ol2, clim1, clim2, llim1, llim2
+pointer buf
+real pmedian, median, dif
+
+int ir_overlap()
+pointer imgs2r()
+real amedr()
+
+begin
+ # Initialize the intensity subraster.
+ call ir_vecinit (deltai, IR_NXSUB(ir) * IR_NYSUB(ir), ranges)
+
+ if (IR_ORDER(ir) == IR_ROW)
+ nmod = IR_NXSUB(ir)
+ else
+ nmod = IR_NYSUB(ir)
+
+ # Loop over the subrasters to be matched.
+ for (num = 1; num <= IR_NXSUB(ir) * IR_NYSUB(ir); num = num + 1) {
+
+ if (num == 1) {
+
+ # Get the position and shift for the first subraster.
+ pideltax = nint (deltax[num])
+ pideltay = nint (deltay[num])
+ pc1 = ic1[num]
+ pc2 = ic2[num]
+ pl1 = il1[num]
+ pl2 = il2[num]
+ num = num + 1
+ dif = 0.0
+ turn_corner = NO
+
+ } else if ((IR_RASTER(ir)) == NO && (mod (num, nmod) == 1)) {
+
+ # Get the position and shift for the first subraster.
+ pideltax = nint (deltax[num-nmod])
+ pideltay = nint (deltay[num-nmod])
+ pc1 = ic1[num-nmod]
+ pc2 = ic2[num-nmod]
+ pl1 = il1[num-nmod]
+ pl2 = il2[num-nmod]
+ dif = -deltai[num-nmod]
+ turn_corner = YES
+
+ } else {
+
+ # Reset the coordinates of the previous subraster.
+ pc1 = c1
+ pc2 = c2
+ pl1 = l1
+ pl2 = l2
+ pideltax = ideltax
+ pideltay = ideltay
+ turn_corner = NO
+ }
+
+ # Get the positions and shifts of the next subraster.
+ ideltax = nint (deltax[num])
+ ideltay = nint (deltay[num])
+ c1 = ic1[num]
+ c2 = ic2[num]
+ l1 = il1[num]
+ l2 = il2[num]
+
+ # Compute the overlap region.
+ if (ir_overlap (pc1 + pideltax, pc2 + pideltax, pl1 + pideltay,
+ pl2 + pideltay, c1 + ideltax, c2 + ideltax, l1 + ideltay,
+ l2 + ideltay, oc1, oc2, ol1, ol2) == YES) {
+
+ clim1 = max (pc1, min (oc1 - pideltax, pc2))
+ clim2 = min (pc2, max (oc2 - pideltax, pc1))
+ llim1 = max (pl1, min (ol1 - pideltay, pl2))
+ llim2 = min (pl2, max (ol2 - pideltay, pl1))
+ buf = imgs2r (im, clim1, clim2, llim1, llim2)
+ pmedian = amedr (Memr[buf], (clim2 - clim1 + 1) * (llim2 -
+ llim1 + 1))
+
+ clim1 = max (c1, min (oc1 - ideltax, c2))
+ clim2 = min (c2, max (oc2 - ideltax, c1))
+ llim1 = max (l1, min (ol1 - ideltay, l2))
+ llim2 = min (l2, max (ol2 - ideltay, l1))
+ buf = imgs2r (im, clim1, clim2, llim1, llim2)
+ median = amedr (Memr[buf], (clim2 - clim1 + 1) * (llim2 -
+ llim1 + 1))
+
+ dif = dif + median - pmedian
+ if (turn_corner == YES) {
+ if (! IS_INDEFR (deltai[num]))
+ deltai[num] = deltai[num-nmod] - median + pmedian
+ } else {
+ if (! IS_INDEFR (deltai[num]))
+ deltai[num] = deltai[num] - dif
+ }
+ }
+
+ }
+end
diff --git a/noao/nproto/ir/irmatch2d.x b/noao/nproto/ir/irmatch2d.x
new file mode 100644
index 00000000..6e1bcd50
--- /dev/null
+++ b/noao/nproto/ir/irmatch2d.x
@@ -0,0 +1,276 @@
+include <imhdr.h>
+include "iralign.h"
+
+# IR_M2MATCH -- Compute the intensity matching parameters.
+
+procedure ir_m2match (ir, im, ranges, ic1, ic2, il1, il2, deltax, deltay,
+ deltai)
+
+pointer ir # pointer to the ir structure
+pointer im # pointer to the input image
+int ranges[ARB] # ranges of data to align
+int ic1[ARB] # array of input begin columns
+int ic2[ARB] # array of input end columns
+int il1[ARB] # array of input begin lines
+int il2[ARB] # array of input end lines
+real deltax[ARB] # array of x shifts
+real deltay[ARB] # array of y shifts
+real deltai[ARB] # array of i shifts
+
+begin
+ # Initialize the intensity subraster.
+ call ir_vecinit (deltai, IR_NXSUB(ir) * IR_NYSUB(ir), ranges)
+ if (ranges[1] == NULL)
+ return
+
+ # Match the intensities in the direction of observation.
+ call ir_omatch (ir, im, ic1, ic2, il1, il2, deltax, deltay, deltai)
+
+ # Match the intensities in the other direction.
+ call ir_nmatch (ir, im, ic1, ic2, il1, il2, deltax, deltay, deltai)
+end
+
+
+# IR_OMATCH -- Procedure to match images in the direction of observation
+# direction.
+
+procedure ir_omatch (ir, im, ic1, ic2, il1, il2, deltax, deltay, deltai)
+
+pointer ir # pointer to the ir structure
+pointer im # pointer to the input image
+int ic1[ARB] # beginning column limits
+int ic2[ARB] # ending column limits
+int il1[ARB] # beginning line limits
+int il2[ARB] # ending line limits
+real deltax[ARB] # array of x shifts
+real deltay[ARB] # array of y shifts
+real deltai[ARB] # array of intensity shifts
+
+int num, nimages, nrasters
+int pc1, pc2, pl1, pl2, c1, c2, l1, l2
+int pideltax, pideltay, ideltax, ideltay
+int oc1, oc2, ol1, ol2, clim1, clim2, llim1, llim2
+pointer buf
+real pmedian, median, dif
+
+int ir_overlap()
+pointer imgs2r()
+real amedr()
+
+begin
+ # Compute the do loop parameters.
+ nimages = IR_NXSUB(ir) * IR_NYSUB(ir)
+ if (IR_ORDER(ir) == IR_ROW)
+ nrasters = IR_NXSUB(ir)
+ else
+ nrasters = IR_NYSUB(ir)
+
+ # Loop over the subrasters to be matched.
+ for (num = 1; num <= nimages; num = num + 1) {
+
+ if (mod (num, nrasters) == 1) {
+
+ # Get the position and shift for the first subraster in
+ # the column.
+ pideltax = nint (deltax[num])
+ pideltay = nint (deltay[num])
+ pc1 = ic1[num]
+ pc2 = ic2[num]
+ pl1 = il1[num]
+ pl2 = il2[num]
+ num = num + 1
+ dif = 0.0
+
+ # Get the the position and shift for the next subraster in
+ # the column.to be
+ ideltax = nint (deltax[num])
+ ideltay = nint (deltay[num])
+ c1 = ic1[num]
+ c2 = ic2[num]
+ l1 = il1[num]
+ l2 = il2[num]
+
+ } else {
+
+ # Reset the coordinates of the previous subraster.
+ pc1 = c1
+ pc2 = c2
+ pl1 = l1
+ pl2 = l2
+ pideltax = ideltax
+ pideltay = ideltay
+
+ # Get the positions and shifts of the next subraster.
+ ideltax = nint (deltax[num])
+ ideltay = nint (deltay[num])
+ c1 = ic1[num]
+ c2 = ic2[num]
+ l1 = il1[num]
+ l2 = il2[num]
+
+ }
+
+ # Compute the overlap region.
+ if (ir_overlap (pc1 + pideltax, pc2 + pideltax, pl1 + pideltay,
+ pl2 + pideltay, c1 + ideltax, c2 + ideltax, l1 + ideltay,
+ l2 + ideltay, oc1, oc2, ol1, ol2) == YES) {
+
+ clim1 = max (pc1, min (oc1 - pideltax, pc2))
+ clim2 = min (pc2, max (oc2 - pideltax, pc1))
+ llim1 = max (pl1, min (ol1 - pideltay, pl2))
+ llim2 = min (pl2, max (ol2 - pideltay, pl1))
+ buf = imgs2r (im, clim1, clim2, llim1, llim2)
+ pmedian = amedr (Memr[buf], (clim2 - clim1 + 1) * (llim2 -
+ llim1 + 1))
+
+ clim1 = max (c1, min (oc1 - ideltax, c2))
+ clim2 = min (c2, max (oc2 - ideltax, c1))
+ llim1 = max (l1, min (ol1 - ideltay, l2))
+ llim2 = min (l2, max (ol2 - ideltay, l1))
+ buf = imgs2r (im, clim1, clim2, llim1, llim2)
+ median = amedr (Memr[buf], (clim2 - clim1 + 1) * (llim2 -
+ llim1 + 1))
+
+ dif = dif + median - pmedian
+ if (! IS_INDEFR (deltai[num]))
+ deltai[num] = deltai[num] - dif
+ }
+ }
+end
+
+
+# IR_NMATCH -- Procedure to match images in the other direction.
+
+procedure ir_nmatch (ir, im, ic1, ic2, il1, il2, deltax, deltay, deltai)
+
+pointer ir # pointer to the ir structure
+pointer im # pointer to the input image
+int ic1[ARB] # array of beginning columns
+int ic2[ARB] # array of ending columns
+int il1[ARB] # array of beginning lines
+int il2[ARB] # array of ending lines
+real deltax[ARB] # array of x shifts
+real deltay[ARB] # array of y shifts
+real deltai[ARB] # array of intensity shifts
+
+int num, nrasters, fac, nimages, count
+int pc1, pc2, pl1, pl2, c1, c2, l1, l2
+int pideltax, pideltay, ideltax, ideltay
+int oc1, oc2, ol1, ol2, clim1, clim2, llim1, llim2
+pointer buf
+real pmedian, median, pdif, dif, tdif
+
+int ir_overlap()
+pointer imgs2r()
+real amedr()
+
+begin
+ # Compute the do loop parameters.
+ nimages = IR_NXSUB(ir) * IR_NYSUB(ir)
+ if (IR_ORDER(ir) == IR_ROW)
+ nrasters = IR_NXSUB(ir)
+ else
+ nrasters = IR_NYSUB(ir)
+ fac = 2 * nrasters
+
+ # Loop over the subrasters to be matched.
+ num = 1
+ count = 1
+ repeat {
+
+ # Get the position and shift for the first subraster.
+ if (num <= nrasters) {
+
+ pideltax = nint (deltax[num])
+ pideltay = nint (deltay[num])
+ pc1 = ic1[num]
+ pc2 = ic2[num]
+ pl1 = il1[num]
+ pl2 = il2[num]
+ if (IS_INDEFR(deltai[num]))
+ pdif = 0.0
+ else
+ pdif = deltai[num]
+ tdif = 0.0
+ if (IR_RASTER(ir) == YES) {
+ num = fac - num + 1
+ fac = fac + fac
+ } else
+ num = num + nrasters
+
+ # Get the the position and shift for the next.
+ ideltax = nint (deltax[num])
+ ideltay = nint (deltay[num])
+ c1 = ic1[num]
+ c2 = ic2[num]
+ l1 = il1[num]
+ l2 = il2[num]
+ if (IS_INDEFR(deltai[num]))
+ dif = 0.0
+ else
+ dif = deltai[num]
+
+ } else {
+
+ # Reset the coordinates of the previous subraster.
+ pc1 = c1
+ pc2 = c2
+ pl1 = l1
+ pl2 = l2
+ pideltax = ideltax
+ pideltay = ideltay
+ pdif = dif
+
+ # Get the positions and shifts of the subraster to be adjusted.
+ ideltax = nint (deltax[num])
+ ideltay = nint (deltay[num])
+ c1 = ic1[num]
+ c2 = ic2[num]
+ l1 = il1[num]
+ l2 = il2[num]
+ if (IS_INDEFR(deltai[num]))
+ dif = 0.0
+ else
+ dif = deltai[num]
+
+ }
+
+ # Compute the overlap region.
+ if (ir_overlap (pc1 + pideltax, pc2 + pideltax, pl1 + pideltay,
+ pl2 + pideltay, c1 + ideltax, c2 + ideltax, l1 + ideltay,
+ l2 + ideltay, oc1, oc2, ol1, ol2) == YES) {
+
+ clim1 = max (pc1, oc1 - pideltax)
+ clim2 = min (pc2, oc2 - pideltax)
+ llim1 = max (pl1, ol1 - pideltay)
+ llim2 = min (pl2, ol2 - pideltay)
+ buf = imgs2r (im, clim1, clim2, llim1, llim2)
+ pmedian = amedr (Memr[buf], (clim2 - clim1 + 1) * (llim2 -
+ llim1 + 1))
+
+ clim1 = max (c1, oc1 - ideltax)
+ clim2 = min (c2, oc2 - ideltax)
+ llim1 = max (l1, ol1 - ideltay)
+ llim2 = min (l2, ol2 - ideltay)
+ buf = imgs2r (im, clim1, clim2, llim1, llim2)
+ median = amedr (Memr[buf], (clim2 - clim1 + 1) * (llim2 -
+ llim1 + 1))
+
+ tdif = tdif + median + dif - pmedian - pdif
+ if (! IS_INDEFR (deltai[num]))
+ deltai[num] = deltai[num] - tdif
+ }
+
+ if (IR_RASTER(ir) == YES) {
+ num = fac - num + 1
+ fac = fac + fac
+ } else
+ num = num + nrasters
+ if (num > nimages) {
+ count = count + 1
+ num = count
+ fac = 2 * nrasters
+ }
+
+ } until (count > nrasters)
+end
diff --git a/noao/nproto/ir/irmedr.x b/noao/nproto/ir/irmedr.x
new file mode 100644
index 00000000..436a3673
--- /dev/null
+++ b/noao/nproto/ir/irmedr.x
@@ -0,0 +1,35 @@
+# IRMEDR -- Procedure to compute the median of an array in which some
+# elements are undefined.
+
+real procedure irmedr (a, aindex, npts)
+
+real a[ARB] # input array
+int aindex[ARB] # definition array
+int npts # number of points
+
+int i, n
+pointer sp, b
+real med
+real asokr()
+
+begin
+ call smark (sp)
+ call salloc (b, npts, TY_REAL)
+
+ n = 0
+ do i = 1, npts {
+ if (aindex[i] > 0) {
+ Memr[b+n] = a[i]
+ n = n + 1
+ }
+ }
+
+ if (n == 0)
+ med = INDEFR
+ else
+ med = asokr (Memr[b], n, (n + 1) / 2)
+
+ call sfree (sp)
+
+ return (med)
+end
diff --git a/noao/nproto/ir/iroverlap.x b/noao/nproto/ir/iroverlap.x
new file mode 100644
index 00000000..822e6e99
--- /dev/null
+++ b/noao/nproto/ir/iroverlap.x
@@ -0,0 +1,40 @@
+# IR_OVERLAP -- Procedure to compute the overlap between two rectangles.
+
+int procedure ir_overlap (pc1out, pc2out, pl1out, pl2out, c1out, c2out,
+ l1out, l2out, oc1out, oc2out, ol1out, ol2out)
+
+int pc1out, pc2out # previous subraster column limits
+int pl1out, pl2out # previous subraster line limits
+int c1out, c2out # current subraster column limits
+int l1out, l2out # current subraster line limits
+int oc1out, oc2out # overlap column limits
+int ol1out, ol2out # overlap line limits
+
+begin
+ # Check for the case where no intersection is present.
+ if (c1out > pc2out || c2out < pc1out || l1out > pl2out ||
+ l2out < pl1out)
+ return (NO)
+
+ # Compute the column overlap limits.
+ if (pc1out <= c1out)
+ oc1out = c1out
+ else
+ oc1out = pc1out
+ if (pc2out <= c2out)
+ oc2out = pc2out
+ else
+ oc2out = c2out
+
+ # Compute the line overlap limits.
+ if (pl1out <= l1out)
+ ol1out = l1out
+ else
+ ol1out = pl1out
+ if (pl2out <= l2out)
+ ol2out = pl2out
+ else
+ ol2out = l2out
+
+ return (YES)
+end
diff --git a/noao/nproto/ir/irqsort.x b/noao/nproto/ir/irqsort.x
new file mode 100644
index 00000000..3c8b710c
--- /dev/null
+++ b/noao/nproto/ir/irqsort.x
@@ -0,0 +1,215 @@
+define LOGPTR 20 # log2(maxpts) (1e6)
+
+# IR_QSORT -- Vector Quicksort. In this version the index array is
+# sorted.
+
+procedure ir_qsortr (data, a, b, npix)
+
+real data[ARB] # data array
+int a[ARB], b[ARB] # index array
+int npix # number of pixels
+
+int i, j, lv[LOGPTR], p, uv[LOGPTR], temp
+real pivot
+
+begin
+ # Initialize the indices for an inplace sort.
+ do i = 1, npix
+ a[i] = i
+ call amovi (a, b, npix)
+
+ p = 1
+ lv[1] = 1
+ uv[1] = npix
+ while (p > 0) {
+
+ # If only one elem in subset pop stack otherwise pivot line.
+ if (lv[p] >= uv[p])
+ p = p - 1
+ else {
+ i = lv[p] - 1
+ j = uv[p]
+ pivot = data[b[j]]
+
+ while (i < j) {
+ for (i=i+1; data[b[i]] < pivot; i=i+1)
+ ;
+ for (j=j-1; j > i; j=j-1)
+ if (data[b[j]] <= pivot)
+ break
+ if (i < j) { # out of order pair
+ temp = b[j] # interchange elements
+ b[j] = b[i]
+ b[i] = temp
+ }
+ }
+
+ j = uv[p] # move pivot to position i
+ temp = b[j] # interchange elements
+ b[j] = b[i]
+ b[i] = temp
+
+ if (i-lv[p] < uv[p] - i) { # stack so shorter done first
+ lv[p+1] = lv[p]
+ uv[p+1] = i - 1
+ lv[p] = i + 1
+ } else {
+ lv[p+1] = i + 1
+ uv[p+1] = uv[p]
+ uv[p] = i - 1
+ }
+
+ p = p + 1 # push onto stack
+ }
+ }
+end
+
+
+# IR_QSORT -- Vector Quicksort. In this version the index array is
+# sorted.
+
+procedure ir_qsorti (data, a, b, npix)
+
+int data[ARB] # data array
+int a[ARB], b[ARB] # index array
+int npix # number of pixels
+
+int i, j, lv[LOGPTR], p, uv[LOGPTR], temp
+int pivot
+
+begin
+ # Initialize the indices for an inplace sort.
+ do i = 1, npix
+ a[i] = i
+ call amovi (a, b, npix)
+
+ p = 1
+ lv[1] = 1
+ uv[1] = npix
+ while (p > 0) {
+
+ # If only one elem in subset pop stack otherwise pivot line.
+ if (lv[p] >= uv[p])
+ p = p - 1
+ else {
+ i = lv[p] - 1
+ j = uv[p]
+ pivot = data[b[j]]
+
+ while (i < j) {
+ for (i=i+1; data[b[i]] < pivot; i=i+1)
+ ;
+ for (j=j-1; j > i; j=j-1)
+ if (data[b[j]] <= pivot)
+ break
+ if (i < j) { # out of order pair
+ temp = b[j] # interchange elements
+ b[j] = b[i]
+ b[i] = temp
+ }
+ }
+
+ j = uv[p] # move pivot to position i
+ temp = b[j] # interchange elements
+ b[j] = b[i]
+ b[i] = temp
+
+ if (i-lv[p] < uv[p] - i) { # stack so shorter done first
+ lv[p+1] = lv[p]
+ uv[p+1] = i - 1
+ lv[p] = i + 1
+ } else {
+ lv[p+1] = i + 1
+ uv[p+1] = uv[p]
+ uv[p] = i - 1
+ }
+
+ p = p + 1 # push onto stack
+ }
+ }
+end
+
+
+# IR_QSORT -- Vector Quicksort. In this version the index array is
+# sorted.
+
+procedure ir_qsortb (data, a, b, npix)
+
+bool data[ARB] # data array
+int a[ARB], b[ARB] # index array
+int npix # number of pixels
+
+int i, j, lv[LOGPTR], p, uv[LOGPTR], temp
+bool pivot
+int ir_compareb()
+
+begin
+ # Initialize the indices for an inplace sort.
+ do i = 1, npix
+ a[i] = i
+ call amovi (a, b, npix)
+
+ p = 1
+ lv[1] = 1
+ uv[1] = npix
+ while (p > 0) {
+
+ # If only one elem in subset pop stack otherwise pivot line.
+ if (lv[p] >= uv[p])
+ p = p - 1
+ else {
+ i = lv[p] - 1
+ j = uv[p]
+ pivot = data[b[j]]
+
+ while (i < j) {
+ #for (i=i+1; data[b[i]] != pivot; i=i+1)
+ for (i=i+1; ir_compareb (data[b[i]], pivot) < 0; i=i+1)
+ ;
+ for (j=j-1; j > i; j=j-1)
+ #if (data[b[j]] != pivot)
+ if (ir_compareb (data[b[j]], pivot) <= 0)
+ break
+ if (i < j) { # out of order pair
+ temp = b[j] # interchange elements
+ b[j] = b[i]
+ b[i] = temp
+ }
+ }
+
+ j = uv[p] # move pivot to position i
+ temp = b[j] # interchange elements
+ b[j] = b[i]
+ b[i] = temp
+
+ if (i-lv[p] < uv[p] - i) { # stack so shorter done first
+ lv[p+1] = lv[p]
+ uv[p+1] = i - 1
+ lv[p] = i + 1
+ } else {
+ lv[p+1] = i + 1
+ uv[p+1] = uv[p]
+ uv[p] = i - 1
+ }
+
+ p = p + 1 # push onto stack
+ }
+ }
+end
+
+
+# IR_COMPAREB -- Compare to booleans for the sort routine.
+
+int procedure ir_compareb (a, b)
+
+bool a # first boolean
+bool b # second boolean
+
+begin
+ if (! a && b)
+ return (-1)
+ else if (a && ! b)
+ return (1)
+ else
+ return (0)
+end
diff --git a/noao/nproto/ir/irtools.x b/noao/nproto/ir/irtools.x
new file mode 100644
index 00000000..50367440
--- /dev/null
+++ b/noao/nproto/ir/irtools.x
@@ -0,0 +1,147 @@
+include <imhdr.h>
+include "iralign.h"
+
+# IR_INIT -- Initialize the ir structure
+
+procedure ir_init (ir)
+
+pointer ir # pointer to the ir strucuture
+
+begin
+ call malloc (ir, LEN_IRSTRUCT, TY_STRUCT)
+
+ IR_IC1(ir) = NULL
+ IR_IC2(ir) = NULL
+ IR_IL1(ir) = NULL
+ IR_IL2(ir) = NULL
+ IR_OC1(ir) = NULL
+ IR_OC2(ir) = NULL
+ IR_OL1(ir) = NULL
+ IR_OL2(ir) = NULL
+ IR_DELTAX(ir) = NULL
+ IR_DELTAY(ir) = NULL
+ IR_DELTAI(ir) = NULL
+ IR_XRSHIFTS(ir) = NULL
+ IR_YRSHIFTS(ir) = NULL
+ IR_NRSHIFTS(ir) = NULL
+ IR_XCSHIFTS(ir) = NULL
+ IR_YCSHIFTS(ir) = NULL
+ IR_NCSHIFTS(ir) = NULL
+end
+
+
+# IR_PARAMS -- Get the ir structure parameters
+
+procedure ir_params (ir, im, outim)
+
+pointer ir # pointer to the ir strucuture
+pointer im # pointer to the input image
+pointer outim # pointer to the output image
+
+int nimcols, nimlines
+real rval
+int clgeti()
+real clgetr()
+
+begin
+ IR_NXRSUB(ir) = clgeti ("nxrsub")
+ if (IS_INDEFI(IR_NXRSUB(ir)) || IR_NXRSUB(ir) < 1 || IR_NXRSUB(ir) >
+ IR_NXSUB(ir))
+ IR_NXRSUB(ir) = (IR_NXSUB(ir) + 1) / 2
+ IR_NYRSUB(ir) = clgeti ("nyrsub")
+ if (IS_INDEFI(IR_NYRSUB(ir)) || IR_NYRSUB(ir) < 1 || IR_NYRSUB(ir) >
+ IR_NYSUB(ir))
+ IR_NYRSUB(ir) = (IR_NYSUB(ir) + 1) / 2
+
+ IR_XREF(ir) = clgeti ("xref")
+ IR_YREF(ir) = clgeti ("yref")
+
+ nimcols = clgeti ("nimcols")
+ if (! IS_INDEFI(nimcols) && nimcols > 0 && nimcols >= IM_LEN(im,1))
+ IM_LEN(outim,1) = nimcols
+ nimlines = clgeti ("nimlines")
+ if (! IS_INDEFI(nimlines) && nimlines > 0 && nimlines >= IM_LEN(im,2))
+ IM_LEN(outim,2) = nimlines
+
+ rval = clgetr ("oval")
+ if (! IS_INDEFR(rval))
+ IR_OVAL(ir) = rval
+end
+
+
+# IR_ARRAYS -- Setup the ir structure arrays.
+
+procedure ir_arrays (ir, nimages)
+
+pointer ir # pointer to the ir strucuture
+int nimages # number of images to be mosaiced
+
+begin
+ call malloc (IR_IC1(ir), nimages, TY_INT)
+ call malloc (IR_IC2(ir), nimages, TY_INT)
+ call malloc (IR_IL1(ir), nimages, TY_INT)
+ call malloc (IR_IL2(ir), nimages, TY_INT)
+ call malloc (IR_OC1(ir), nimages, TY_INT)
+ call malloc (IR_OC2(ir), nimages, TY_INT)
+ call malloc (IR_OL1(ir), nimages, TY_INT)
+ call malloc (IR_OL2(ir), nimages, TY_INT)
+ call malloc (IR_DELTAX(ir), nimages, TY_REAL)
+ call malloc (IR_DELTAY(ir), nimages, TY_REAL)
+ call malloc (IR_DELTAI(ir), nimages, TY_REAL)
+
+ call malloc (IR_XRSHIFTS(ir), nimages, TY_REAL)
+ call malloc (IR_YRSHIFTS(ir), nimages, TY_REAL)
+ call malloc (IR_NRSHIFTS(ir), nimages, TY_INT)
+ call malloc (IR_XCSHIFTS(ir), nimages, TY_REAL)
+ call malloc (IR_YCSHIFTS(ir), nimages, TY_REAL)
+ call malloc (IR_NCSHIFTS(ir), nimages, TY_INT)
+end
+
+
+# IR_FREE -- Free the ir strucuture.
+
+procedure ir_free (ir)
+
+pointer ir # pointer to the ir strucuture
+
+begin
+ if (IR_IC1(ir) != NULL)
+ call mfree (IR_IC1(ir), TY_INT)
+ if (IR_IC2(ir) != NULL)
+ call mfree (IR_IC2(ir), TY_INT)
+ if (IR_IL1(ir) != NULL)
+ call mfree (IR_IL1(ir), TY_INT)
+ if (IR_IL2(ir) != NULL)
+ call mfree (IR_IL2(ir), TY_INT)
+ if (IR_OC1(ir) != NULL)
+ call mfree (IR_OC1(ir), TY_INT)
+ if (IR_OC2(ir) != NULL)
+ call mfree (IR_OC2(ir), TY_INT)
+ if (IR_OL1(ir) != NULL)
+ call mfree (IR_OL1(ir), TY_INT)
+ if (IR_OL2(ir) != NULL)
+ call mfree (IR_OL2(ir), TY_INT)
+
+ if (IR_DELTAX(ir) != NULL)
+ call mfree (IR_DELTAX(ir), TY_REAL)
+ if (IR_DELTAY(ir) != NULL)
+ call mfree (IR_DELTAY(ir), TY_REAL)
+ if (IR_DELTAI(ir) != NULL)
+ call mfree (IR_DELTAI(ir), TY_REAL)
+
+ if (IR_XRSHIFTS(ir) != NULL)
+ call mfree (IR_XRSHIFTS(ir), TY_REAL)
+ if (IR_YRSHIFTS(ir) != NULL)
+ call mfree (IR_YRSHIFTS(ir), TY_REAL)
+ if (IR_NRSHIFTS(ir) != NULL)
+ call mfree (IR_NRSHIFTS(ir), TY_INT)
+ if (IR_XCSHIFTS(ir) != NULL)
+ call mfree (IR_XCSHIFTS(ir), TY_REAL)
+ if (IR_YCSHIFTS(ir) != NULL)
+ call mfree (IR_YCSHIFTS(ir), TY_REAL)
+ if (IR_NCSHIFTS(ir) != NULL)
+ call mfree (IR_NCSHIFTS(ir), TY_INT)
+
+ if (ir != NULL)
+ call mfree (ir, TY_STRUCT)
+end
diff --git a/noao/nproto/ir/mkpkg b/noao/nproto/ir/mkpkg
new file mode 100644
index 00000000..7297fa37
--- /dev/null
+++ b/noao/nproto/ir/mkpkg
@@ -0,0 +1,24 @@
+$checkout libpkg.a ../
+$update libpkg.a
+$checkin libpkg.a ../
+$exit
+
+libpkg.a:
+ iralign.x "iralign.h" <imhdr.h>
+ irdbio.x "iralign.h"
+ iriinit.x
+ irindices.x "iralign.h"
+ irimisec.x <ctype.h>
+ irimzero.x
+ irmatch1d.x "iralign.h" <imhdr.h> <pkg/dttext.h>
+ irmatch2d.x "iralign.h" <imhdr.h>
+ irmedr.x
+ iroverlap.x
+ irqsort.x
+ irlinks.x "iralign.h"
+ irtools.x "iralign.h" <imhdr.h>
+ t_iralign.x "iralign.h" <imhdr.h> <fset.h>
+ t_irmatch1d.x "iralign.h" <imhdr.h> <fset.h>
+ t_irmatch2d.x "iralign.h" <imhdr.h> <fset.h>
+ t_irmosaic.x "iralign.h" <imhdr.h> <fset.h>
+ ;
diff --git a/noao/nproto/ir/t_iralign.x b/noao/nproto/ir/t_iralign.x
new file mode 100644
index 00000000..03f97b32
--- /dev/null
+++ b/noao/nproto/ir/t_iralign.x
@@ -0,0 +1,134 @@
+include <imhdr.h>
+include <fset.h>
+include "iralign.h"
+
+# T_IRALIGN -- Align the individual subraster elements in the input image.
+# In order to run this program the user should have created the output image
+# and the database file with the IRMOSAIC task. In addition the user should
+# supply a coordinate list consisting of pairs of coordinates of identical
+# objects or features in two adjacent subrasters.
+
+procedure t_iralign ()
+
+int cl, nimages, interp, align, verbose
+pointer ir, sp, inimage, outimage, database, coords, trimlimits, str
+pointer im, outim, dt
+
+bool clgetb()
+int open(), clgwrd(), btoi()
+int ir_links(), ir_clinks(), ir_flinks()
+pointer immap(), dtmap()
+real clgetr()
+
+begin
+ # Allocate sapce for the ir strucuture.
+ call ir_init (ir)
+
+ # Set the standard output to flush on a new line.
+ call fseti (STDOUT, F_FLUSHNL, YES)
+
+ # Allocate temporary working space.
+ call smark (sp)
+ call salloc (inimage, SZ_FNAME, TY_CHAR)
+ call salloc (outimage, SZ_FNAME, TY_CHAR)
+ call salloc (coords, SZ_FNAME, TY_CHAR)
+ call salloc (database, SZ_FNAME, TY_CHAR)
+ call salloc (trimlimits, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Get the input and output images and the coordinate list.
+ call clgstr ("input", Memc[inimage], SZ_FNAME)
+ call clgstr ("output", Memc[outimage], SZ_FNAME)
+ call clgstr ("database", Memc[database], SZ_FNAME)
+ call clgstr ("coords", Memc[coords], SZ_FNAME)
+ align = clgwrd ("alignment", Memc[str], SZ_LINE, ",coords,shifts,file,")
+ call clgstr ("trimlimits", Memc[trimlimits], SZ_FNAME)
+
+ # Open the images and files.
+ im = immap (Memc[inimage], READ_ONLY, 0)
+ outim = immap (Memc[outimage], NEW_COPY, im)
+ dt = dtmap (Memc[database], READ_ONLY)
+
+ # Get the data base parameters.
+ call ir_dtrparams (dt, Memc[inimage], ir)
+
+ # Get the rest of the parameters.
+ call ir_params (ir, im, outim)
+ interp = clgwrd ("interpolant", Memc[str], SZ_LINE,
+ ",nearest,linear,poly3,poly5,spline3,")
+ verbose = btoi (clgetb ("verbose"))
+
+ # Allocate array space.
+ nimages = IR_NXSUB(ir) * IR_NYSUB(ir)
+ call ir_arrays (ir, nimages)
+
+ # Compute the shifts for each subraster.
+ switch (align) {
+ case IR_COORDS:
+ cl = open (Memc[coords], READ_ONLY, TEXT_FILE)
+ if (ir_links (cl, Memr[IR_XRSHIFTS(ir)], Memr[IR_YRSHIFTS(ir)],
+ Memr[IR_XCSHIFTS(ir)], Memr[IR_YCSHIFTS(ir)],
+ Memi[IR_NRSHIFTS(ir)], Memi[IR_NCSHIFTS(ir)],
+ IR_NCOLS(ir), IR_NROWS(ir), IR_NXRSUB(ir), IR_NYRSUB(ir),
+ IR_NXSUB(ir), IR_NYSUB(ir), IR_NXOVERLAP(ir), IR_NYOVERLAP(ir),
+ IR_ORDER(ir)) > 0) {
+ call ir_shifts (ir, im, outim, Memr[IR_XRSHIFTS(ir)],
+ Memr[IR_YRSHIFTS(ir)], Memr[IR_XCSHIFTS(ir)],
+ Memr[IR_YCSHIFTS(ir)], Memi[IR_IC1(ir)],
+ Memi[IR_IC2(ir)], Memi[IR_IL1(ir)], Memi[IR_IL2(ir)],
+ Memi[IR_OC1(ir)], Memi[IR_OC2(ir)], Memi[IR_OL1(ir)],
+ Memi[IR_OL2(ir)], Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)])
+ } else
+ call error (0, "There are no legal shifts in the coords file.")
+ call close (cl)
+
+ case IR_SHIFTS:
+ if (ir_clinks (Memr[IR_XRSHIFTS(ir)], Memr[IR_YRSHIFTS(ir)],
+ Memr[IR_XCSHIFTS(ir)], Memr[IR_YCSHIFTS(Ir)], IR_NXRSUB(ir),
+ IR_NYRSUB(ir), IR_NXSUB(ir), IR_NYSUB(ir), clgetr ("xshift"),
+ clgetr ("yshift")) > 0) {
+ call ir_shifts (ir, im, outim, Memr[IR_XRSHIFTS(ir)],
+ Memr[IR_YRSHIFTS(ir)], Memr[IR_XCSHIFTS(ir)],
+ Memr[IR_YCSHIFTS(ir)], Memi[IR_IC1(ir)],
+ Memi[IR_IC2(ir)], Memi[IR_IL1(ir)], Memi[IR_IL2(ir)],
+ Memi[IR_OC1(ir)], Memi[IR_OC2(ir)], Memi[IR_OL1(ir)],
+ Memi[IR_OL2(ir)], Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)])
+ } else
+ call error (0, "There are no legal shifts in the coords file.")
+
+ case IR_FILE:
+
+ cl = open (Memc[coords], READ_ONLY, TEXT_FILE)
+ if (ir_flinks (cl, Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)],
+ Memr[IR_DELTAI(ir)], nimages) >= nimages) {
+ call ir_fshifts (ir, im, outim, Memr[IR_DELTAX(ir)],
+ Memr[IR_DELTAY(ir)], Memi[IR_IC1(ir)], Memi[IR_IC2(ir)],
+ Memi[IR_IL1(ir)], Memi[IR_IL2(ir)], Memi[IR_OC1(ir)],
+ Memi[IR_OC2(ir)], Memi[IR_OL1(ir)], Memi[IR_OL2(ir)])
+ } else
+ call error (0, "There are fewer shifts than subraster.")
+ call close (cl)
+
+ default:
+ call error (0, "T_IRALIGN: Undefined alignment algorithm")
+ }
+
+ # Fill the output image with the undefined value.
+ call ir_imzero (outim, int (IM_LEN(outim,1)), int (IM_LEN(outim,2)),
+ IR_OVAL(ir))
+
+ # Shift all the subrasters.
+ call amovkr (0.0, Memr[IR_DELTAI(ir)], nimages)
+ call ir_subalign (ir, im, outim, Memc[trimlimits], Memi[IR_IC1(ir)],
+ Memi[IR_IC2(ir)], Memi[IR_IL1(ir)], Memi[IR_IL2(ir)],
+ Memi[IR_OC1(ir)], Memi[IR_OC2(ir)], Memi[IR_OL1(ir)],
+ Memi[IR_OL2(ir)], Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)],
+ Memr[IR_DELTAI(ir)], NO, interp, verbose)
+
+ # Close up files
+ call imunmap (im)
+ call imunmap (outim)
+ call dtunmap (dt)
+ call sfree (sp)
+ call ir_free (ir)
+end
diff --git a/noao/nproto/ir/t_irmatch1d.x b/noao/nproto/ir/t_irmatch1d.x
new file mode 100644
index 00000000..05e89e0a
--- /dev/null
+++ b/noao/nproto/ir/t_irmatch1d.x
@@ -0,0 +1,159 @@
+include <imhdr.h>
+include <fset.h>
+include "iralign.h"
+
+# T_IRMATCHD1 -- Align the individual subraster elements in the input image.
+# In order to run this program the user should have created the output image
+# and the database file with the IRMOSAIC task. In addition the user should
+# supply a coordinate list consisting of pairs of coordinates of identical
+# objects or features in two adjacent subrasters.
+
+procedure t_irmatchd1 ()
+
+int cl, interp, align, verbose, nmatch, nimages
+pointer sp, inimage, outimage, database, coords, matchlist, trimlimits, ranges
+pointer str, ir, im, outim, dt
+
+bool clgetb()
+int open(), clgwrd(), btoi(), ir_links(), ir_clinks, ir_flinks()
+int decode_ranges()
+pointer immap(), dtmap()
+real clgetr()
+
+begin
+ # Allocate space for the ir strucuture.
+ call ir_init (ir)
+
+ # Set the standard output to flush on a new line.
+ call fseti (STDOUT, F_FLUSHNL, YES)
+
+ # Allocate temporary working space.
+ call smark (sp)
+ call salloc (inimage, SZ_FNAME, TY_CHAR)
+ call salloc (outimage, SZ_FNAME, TY_CHAR)
+ call salloc (coords, SZ_FNAME, TY_CHAR)
+ call salloc (database, SZ_FNAME, TY_CHAR)
+ call salloc (matchlist, SZ_LINE, TY_CHAR)
+ call salloc (trimlimits, SZ_FNAME, TY_CHAR)
+ call salloc (ranges, 3 * MAX_NRANGES + 1, TY_INT)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Get the input and output images and the coordinate list.
+ call clgstr ("input", Memc[inimage], SZ_FNAME)
+ call clgstr ("output", Memc[outimage], SZ_FNAME)
+ call clgstr ("database", Memc[database], SZ_FNAME)
+ call clgstr ("coords", Memc[coords], SZ_FNAME)
+ align = clgwrd ("alignment", Memc[str], SZ_LINE, ",coords,shifts,file,")
+ call clgstr ("match", Memc[matchlist], SZ_LINE)
+ call clgstr ("trimlimits", Memc[trimlimits], SZ_LINE)
+
+ # Open the images and files.
+ im = immap (Memc[inimage], READ_ONLY, 0)
+ outim = immap (Memc[outimage], NEW_COPY, im)
+ dt = dtmap (Memc[database], READ_ONLY)
+
+ # Get the data base parameters.
+ call ir_dtrparams (dt, Memc[inimage], ir)
+
+ # Get the remaining parameters.
+ call ir_params (ir, im, outim)
+
+ interp = clgwrd ("interpolant", Memc[str], SZ_LINE,
+ ",nearest,linear,poly3,poly5,spline3,")
+ verbose = btoi (clgetb ("verbose"))
+
+ # Decode the list of input images to be intensity matched.
+ nimages = IR_NXSUB(ir) * IR_NYSUB(ir)
+ if (Memc[matchlist] == EOS) {
+ Memi[ranges] = NULL
+ } else if (Memc[matchlist] == '*') {
+ Memi[ranges] = 1
+ Memi[ranges+1] = nimages
+ Memi[ranges+2] = 1
+ Memi[ranges+3] = NULL
+ } else if (decode_ranges (Memc[matchlist], Memi[ranges], MAX_NRANGES,
+ nmatch) == ERR) {
+ call error (0,
+ "Cannot decode list of rasters to be intensity matched.")
+
+ }
+
+ # Allocate space for the ir arrays.
+ call ir_arrays (ir, nimages)
+
+ # Compute the shifts for each subraster.
+ switch (align) {
+ case IR_COORDS:
+ cl = open (Memc[coords], READ_ONLY, TEXT_FILE)
+ if (ir_links (cl, Memr[IR_XRSHIFTS(ir)], Memr[IR_YRSHIFTS(ir)],
+ Memr[IR_XCSHIFTS(ir)], Memr[IR_YCSHIFTS(ir)],
+ Memi[IR_NRSHIFTS(ir)], Memi[IR_NCSHIFTS(ir)],
+ IR_NCOLS(ir), IR_NROWS(ir), IR_NXRSUB(ir), IR_NYRSUB(ir),
+ IR_NXSUB(ir), IR_NYSUB(ir), IR_NXOVERLAP(ir), IR_NYOVERLAP(ir),
+ IR_ORDER(ir)) > 0) {
+ call ir_shifts (ir, im, outim, Memr[IR_XRSHIFTS(ir)],
+ Memr[IR_YRSHIFTS(ir)], Memr[IR_XCSHIFTS(ir)],
+ Memr[IR_YCSHIFTS(ir)], Memi[IR_IC1(ir)], Memi[IR_IC2(ir)],
+ Memi[IR_IL1(ir)], Memi[IR_IL2(ir)], Memi[IR_OC1(ir)],
+ Memi[IR_OC2(ir)], Memi[IR_OL1(ir)], Memi[IR_OL2(ir)],
+ Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)])
+ call ir_m1match (ir, im, Memi[ranges], Memi[IR_IC1(ir)],
+ Memi[IR_IC2(ir)], Memi[IR_IL1(ir)], Memi[IR_IL2(ir)],
+ Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)],
+ Memr[IR_DELTAI(ir)])
+ } else
+ call error (0, "There are no legal shifts in the coords file.")
+ call close (cl)
+
+ case IR_SHIFTS:
+ if (ir_clinks (Memr[IR_XRSHIFTS(ir)], Memr[IR_YRSHIFTS(ir)],
+ Memr[IR_XCSHIFTS(ir)], Memr[IR_YCSHIFTS(ir)], IR_NXRSUB(ir),
+ IR_NYRSUB(ir), IR_NXSUB(ir), IR_NYSUB(ir), clgetr ("xshift"),
+ clgetr ("yshift")) > 0) {
+ call ir_shifts (ir, im, outim, Memr[IR_XRSHIFTS(ir)],
+ Memr[IR_YRSHIFTS(ir)], Memr[IR_XCSHIFTS(ir)],
+ Memr[IR_YCSHIFTS(ir)], Memi[IR_IC1(ir)], Memi[IR_IC2(ir)],
+ Memi[IR_IL1(ir)], Memi[IR_IL2(ir)], Memi[IR_OC1(ir)],
+ Memi[IR_OC2(ir)], Memi[IR_OL1(ir)], Memi[IR_OL2(ir)],
+ Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)])
+ call ir_m1match (ir, im, Memi[ranges], Memi[IR_IC1(ir)],
+ Memi[IR_IC2(ir)], Memi[IR_IL1(ir)], Memi[IR_IL2(ir)],
+ Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)],
+ Memr[IR_DELTAI(ir)])
+ } else
+ call error (0, "There are no legal shifts in the coords file.")
+
+ case IR_FILE:
+ cl = open (Memc[coords], READ_ONLY, TEXT_FILE)
+ if (ir_flinks (cl, Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)],
+ Memr[IR_DELTAI(ir)], nimages) >= nimages) {
+ call ir_fshifts (ir, im, outim, Memr[IR_DELTAX(ir)],
+ Memr[IR_DELTAY(ir)], Memi[IR_IC1(ir)], Memi[IR_IC2(ir)],
+ Memi[IR_IL1(ir)], Memi[IR_IL2(ir)], Memi[IR_OC1(ir)],
+ Memi[IR_OC2(ir)], Memi[IR_OL1(ir)], Memi[IR_OL2(ir)])
+ } else
+ call error (0, "There are fewer shifts than subrasters.")
+ call close (cl)
+
+ default:
+ call error (0, "T_IRALIGN: Undefined alignment algorithm")
+ }
+
+ # Fill the output image with undefined values.
+ call ir_imzero (outim, int (IM_LEN(outim,1)), int (IM_LEN(outim,2)),
+ IR_OVAL(ir))
+
+ # Shift all the subrasters.
+ call ir_subalign (ir, im, outim, Memc[trimlimits], Memi[IR_IC1(ir)],
+ Memi[IR_IC2(ir)], Memi[IR_IL1(ir)], Memi[IR_IL2(ir)],
+ Memi[IR_OC1(ir)], Memi[IR_OC2(ir)], Memi[IR_OL1(ir)],
+ Memi[IR_OL2(ir)], Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)],
+ Memr[IR_DELTAI(ir)], YES, interp, verbose)
+
+ # Close up files
+ call imunmap (im)
+ call imunmap (outim)
+ call dtunmap (dt)
+ call sfree (sp)
+ call ir_free (ir)
+end
diff --git a/noao/nproto/ir/t_irmatch2d.x b/noao/nproto/ir/t_irmatch2d.x
new file mode 100644
index 00000000..38c6feb3
--- /dev/null
+++ b/noao/nproto/ir/t_irmatch2d.x
@@ -0,0 +1,159 @@
+include <imhdr.h>
+include <fset.h>
+include "iralign.h"
+
+# T_IRMATCHD2 -- Align the individual subraster elements in the input image.
+# In order to run this program the user should have created the output image
+# and the database file with the IRMOSAIC task. In addition the user should
+# supply a coordinate list consisting of pairs of coordinates of identical
+# objects or features in two adjacent subrasters.
+
+procedure t_irmatchd2 ()
+
+int cl, interp, align, verbose, nimages, nmatch
+pointer sp, inimage, outimage, database, coords, matchlist, trimlimits, ranges
+pointer str, ir, im, outim, dt
+
+bool clgetb()
+int open(), clgwrd(), btoi(), ir_links(), ir_clinks(), ir_flinks()
+int decode_ranges()
+pointer immap(), dtmap()
+real clgetr()
+
+begin
+ # Allocate space for the ir strucuture.
+ call ir_init (ir)
+
+ # Set the standard output to flush on a new line.
+ call fseti (STDOUT, F_FLUSHNL, YES)
+
+ # Allocate temporary working space.
+ call smark (sp)
+ call salloc (inimage, SZ_FNAME, TY_CHAR)
+ call salloc (outimage, SZ_FNAME, TY_CHAR)
+ call salloc (coords, SZ_FNAME, TY_CHAR)
+ call salloc (database, SZ_FNAME, TY_CHAR)
+ call salloc (matchlist, SZ_LINE, TY_CHAR)
+ call salloc (trimlimits, SZ_FNAME, TY_CHAR)
+ call salloc (ranges, 3 * MAX_NRANGES + 1, TY_INT)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Get the input and output images and the coordinate list.
+ call clgstr ("input", Memc[inimage], SZ_FNAME)
+ call clgstr ("output", Memc[outimage], SZ_FNAME)
+ call clgstr ("database", Memc[database], SZ_FNAME)
+ call clgstr ("coords", Memc[coords], SZ_FNAME)
+ align = clgwrd ("alignment", Memc[str], SZ_LINE, ",coords,shifts,file,")
+ call clgstr ("match", Memc[matchlist], SZ_LINE)
+ call clgstr ("trimlimits", Memc[trimlimits], SZ_FNAME)
+
+ # Open the images and files.
+ im = immap (Memc[inimage], READ_ONLY, 0)
+ outim = immap (Memc[outimage], NEW_COPY, im)
+ dt = dtmap (Memc[database], READ_ONLY)
+
+ # Get the data base parameters.
+ call ir_dtrparams (dt, Memc[inimage], ir)
+
+ call ir_params (ir, im, outim)
+
+ interp = clgwrd ("interpolant", Memc[str], SZ_LINE,
+ ",nearest,linear,poly3,poly5,spline3,")
+ verbose = btoi (clgetb ("verbose"))
+
+ # Decode the list of input images to be intensity matched.
+ nimages = IR_NXSUB(ir) * IR_NYSUB(ir)
+ if (Memc[matchlist] == EOS) {
+ Memi[ranges] = NULL
+ } else if (Memc[matchlist] == '*') {
+ Memi[ranges] = 1
+ Memi[ranges+1] = nimages
+ Memi[ranges+2] = 1
+ Memi[ranges+3] = NULL
+ } else if (decode_ranges (Memc[matchlist], Memi[ranges], MAX_NRANGES,
+ nmatch) == ERR) {
+ call error (0,
+ "Cannot decode list of rasters to be intensity matched.")
+
+ }
+
+ # Allocate working space.
+ call ir_arrays (ir, nimages)
+
+ # Compute the shifts for each subraster.
+ switch (align) {
+ case IR_COORDS:
+ cl = open (Memc[coords], READ_ONLY, TEXT_FILE)
+ if (ir_links (cl, Memr[IR_XRSHIFTS(ir)], Memr[IR_YRSHIFTS(ir)],
+ Memr[IR_XCSHIFTS(ir)], Memr[IR_YCSHIFTS(ir)],
+ Memi[IR_NRSHIFTS(ir)], Memi[IR_NCSHIFTS(ir)],
+ IR_NCOLS(ir), IR_NROWS(ir), IR_NXRSUB(ir), IR_NYRSUB(ir),
+ IR_NXSUB(ir), IR_NYSUB(ir), IR_NXOVERLAP(ir), IR_NYOVERLAP(ir),
+ IR_ORDER(ir)) > 0) {
+ call ir_shifts (ir, im, outim, Memr[IR_XRSHIFTS(ir)],
+ Memr[IR_YRSHIFTS(ir)], Memr[IR_XCSHIFTS(ir)],
+ Memr[IR_YCSHIFTS(ir)], Memi[IR_IC1(ir)], Memi[IR_IC2(ir)],
+ Memi[IR_IL1(ir)], Memi[IR_IL2(ir)], Memi[IR_OC1(ir)],
+ Memi[IR_OC2(ir)], Memi[IR_OL1(ir)], Memi[IR_OL2(ir)],
+ Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)])
+ call ir_m2match (ir, im, Memi[ranges], Memi[IR_IC1(ir)],
+ Memi[IR_IC2(ir)], Memi[IR_IL1(ir)], Memi[IR_IL2(ir)],
+ Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)],
+ Memr[IR_DELTAI(ir)])
+ } else
+ call error (0, "There are no legal shifts in the coords file.")
+ call close (cl)
+
+ case IR_SHIFTS:
+ if (ir_clinks (Memr[IR_XRSHIFTS(ir)], Memr[IR_YRSHIFTS(ir)],
+ Memr[IR_XCSHIFTS(ir)], Memr[IR_YCSHIFTS(ir)], IR_NXRSUB(ir),
+ IR_NYRSUB(ir), IR_NXSUB(ir), IR_NYSUB(ir), clgetr ("xshift"),
+ clgetr ("yshift")) > 0) {
+ call ir_shifts (ir, im, outim, Memr[IR_XRSHIFTS(ir)],
+ Memr[IR_YRSHIFTS(ir)], Memr[IR_XCSHIFTS(ir)],
+ Memr[IR_YCSHIFTS(ir)], Memi[IR_IC1(ir)], Memi[IR_IC2(ir)],
+ Memi[IR_IL1(ir)], Memi[IR_IL2(ir)], Memi[IR_OC1(ir)],
+ Memi[IR_OC2(ir)], Memi[IR_OL1(ir)], Memi[IR_OL2(ir)],
+ Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)])
+ call ir_m2match (ir, im, Memi[ranges], Memi[IR_IC1(ir)],
+ Memi[IR_IC2(ir)], Memi[IR_IL1(ir)], Memi[IR_IL2(ir)],
+ Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)],
+ Memr[IR_DELTAI(ir)])
+ } else
+ call error (0, "There are no legal shifts in the coords file.")
+
+ case IR_FILE:
+ cl = open (Memc[coords], READ_ONLY, TEXT_FILE)
+ if (ir_flinks (cl, Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)],
+ Memr[IR_DELTAI(ir)], nimages) >= nimages) {
+ call ir_fshifts (ir, im, outim, Memr[IR_DELTAX(ir)],
+ Memr[IR_DELTAY(ir)], Memi[IR_IC1(ir)], Memi[IR_IC2(ir)],
+ Memi[IR_IL1(ir)], Memi[IR_IL2(ir)], Memi[IR_OC1(ir)],
+ Memi[IR_OC2(ir)], Memi[IR_OL1(ir)], Memi[IR_OL2(ir)])
+ } else
+ call error (0, "There are fewer shifts than input subrasters.")
+ call close (cl)
+
+ default:
+ call error (0, "T_IRALIGN: Undefined alignment algorithm")
+ }
+
+
+ # Fill the output image with the unknown value.
+ call ir_imzero (outim, int (IM_LEN(outim,1)), int (IM_LEN(outim, 2)),
+ IR_OVAL(ir))
+
+ # Shift and match all the subrasters.
+ call ir_subalign (ir, im, outim, Memc[trimlimits], Memi[IR_IC1(ir)],
+ Memi[IR_IC2(ir)], Memi[IR_IL1(ir)], Memi[IR_IL2(ir)],
+ Memi[IR_OC1(ir)], Memi[IR_OC2(ir)], Memi[IR_OL1(ir)],
+ Memi[IR_OL2(ir)], Memr[IR_DELTAX(ir)], Memr[IR_DELTAY(ir)],
+ Memr[IR_DELTAI(ir)], YES, interp, verbose)
+
+ # Close up files.
+ call imunmap (im)
+ call imunmap (outim)
+ call dtunmap (dt)
+ call sfree (sp)
+ call ir_free (ir)
+end
diff --git a/noao/nproto/ir/t_irmosaic.x b/noao/nproto/ir/t_irmosaic.x
new file mode 100644
index 00000000..86dee342
--- /dev/null
+++ b/noao/nproto/ir/t_irmosaic.x
@@ -0,0 +1,498 @@
+include <imhdr.h>
+include <fset.h>
+include "iralign.h"
+
+
+# T_IRMOSAIC -- Procedure to combine a list of subrasters into a single large
+# image.
+
+procedure t_irmosaic ()
+
+int nimages, nmissing, verbose, subtract
+pointer ir, sp, outimage, database, trimsection, medsection, nullinput, ranges
+pointer str, index, c1, c2, l1, l2, isnull, median, imlist, outim, dt
+
+bool clgetb()
+char clgetc()
+int btoi(), clgwrd(), imtlen(), clgeti(), decode_ranges(), ir_get_imtype()
+pointer imtopenp(), ir_setim(), dtmap()
+real clgetr()
+
+begin
+ call fseti (STDOUT, F_FLUSHNL, YES)
+ call malloc (ir, LEN_IRSTRUCT, TY_STRUCT)
+
+ # Allocate temporary working space.
+ call smark (sp)
+ call salloc (outimage, SZ_FNAME, TY_CHAR)
+ call salloc (database, SZ_FNAME, TY_CHAR)
+ call salloc (trimsection, SZ_FNAME, TY_CHAR)
+ call salloc (medsection, SZ_FNAME, TY_CHAR)
+ call salloc (nullinput, SZ_FNAME, TY_CHAR)
+ call salloc (ranges, 3 * MAX_NRANGES + 1, TY_INT)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Get the image list, output image name and database file name.
+ imlist = imtopenp ("input")
+ call clgstr ("output", Memc[outimage], SZ_FNAME)
+ call clgstr ("database", Memc[database], SZ_FNAME)
+ call clgstr ("trim_section", Memc[trimsection], SZ_FNAME)
+ call clgstr ("null_input", Memc[nullinput], SZ_FNAME)
+ call clgstr ("median_section", Memc[medsection], SZ_FNAME)
+ if (Memc[medsection] == EOS)
+ subtract = NO
+ else
+ subtract = btoi (clgetb ("subtract"))
+ verbose = btoi (clgetb ("verbose"))
+
+ # Get the mosaicing parameters.
+ IR_NXSUB(ir) = clgeti ("nxsub")
+ IR_NYSUB(ir) = clgeti ("nysub")
+ IR_CORNER(ir) = clgwrd ("corner", Memc[str], SZ_FNAME, ",ll,lr,ul,ur,")
+ IR_ORDER(ir) = clgwrd ("direction", Memc[str], SZ_FNAME, ",row,column,")
+ IR_RASTER(ir) = btoi (clgetb ("raster"))
+ IR_NXOVERLAP(ir) = clgeti ("nxoverlap")
+ IR_NYOVERLAP(ir) = clgeti ("nyoverlap")
+ IR_OVAL(ir) = clgetr ("oval")
+
+ # Check that the number of observed and missing images matches
+ # the number of specified subrasters.
+ if (Memc[nullinput] == EOS) {
+ nmissing = 0
+ Memi[ranges] = 0
+ Memi[ranges+1] = 0
+ Memi[ranges+2] = 1
+ Memi[ranges+3] = NULL
+ } else {
+ if (decode_ranges (Memc[nullinput], Memi[ranges], MAX_NRANGES,
+ nmissing) == ERR)
+ call error (0, "Error decoding list of unobserved rasters.")
+ }
+ nimages = imtlen (imlist) + nmissing
+ if (nimages != (IR_NXSUB(ir) * IR_NYSUB(ir)))
+ call error (0,
+ "The number of input images is not equal to nxsub * nysub.")
+
+ # Compute the output image characteristics and open the output image.
+ outim = ir_setim (ir, imlist, Memc[trimsection], Memc[outimage],
+ clgeti ("nimcols"), clgeti ("nimrows"), ir_get_imtype (clgetc (
+ "opixtype")))
+
+ # Open the database file.
+ dt = dtmap (Memc[database], APPEND)
+
+ # Allocate space for and setup the database.
+ call salloc (index, nimages, TY_INT)
+ call salloc (c1, nimages, TY_INT)
+ call salloc (c2, nimages, TY_INT)
+ call salloc (l1, nimages, TY_INT)
+ call salloc (l2, nimages, TY_INT)
+ call salloc (isnull, nimages, TY_INT)
+ call salloc (median, nimages, TY_REAL)
+
+ call ir_setup (ir, imlist, Memi[ranges], Memc[trimsection],
+ Memc[medsection], outim, Memi[index], Memi[c1], Memi[c2],
+ Memi[l1], Memi[l2], Memi[isnull], Memr[median])
+
+ # Write the parameters to the database file.
+ call ir_dtwparams (dt, Memc[outimage], Memc[trimsection],
+ Memc[medsection], ir)
+
+ # Make the output image.
+ call ir_mkmosaic (imlist, Memc[trimsection], outim, Memi[index],
+ Memi[c1], Memi[c2], Memi[l1], Memi[l2], Memi[isnull],
+ Memr[median], IR_NXSUB(ir), IR_NYSUB(ir), IR_OVAL(ir), subtract)
+
+ # Write the database file.
+ call ir_dtwinput (imlist, Memc[trimsection], Memc[outimage], dt,
+ Memi[index], Memi[c1], Memi[c2], Memi[l1], Memi[l2], Memi[isnull],
+ Memr[median], IR_NXSUB(ir) * IR_NYSUB(ir), subtract, verbose)
+
+ # Close up files and free space.
+ call dtunmap (dt)
+ call imunmap (outim)
+ call clpcls (imlist)
+ call sfree (sp)
+ call mfree (ir, TY_STRUCT)
+end
+
+
+define NTYPES 7
+
+# IR_GET_IMTYPE -- Procedure to get the image type.
+
+int procedure ir_get_imtype (c)
+
+char c # character denoting the image type
+
+int i, typecodes[NTYPES]
+int stridx()
+string types "usilrdx"
+data typecodes /TY_USHORT, TY_SHORT, TY_INT, TY_LONG, TY_REAL, TY_DOUBLE,
+ TY_COMPLEX/
+
+begin
+ i = stridx (c, types)
+ if (i == 0)
+ return (ERR)
+ else
+ return (typecodes[i])
+end
+
+
+# IR_SETUP -- Setup the data base parameters for the images.
+
+procedure ir_setup (ir, imlist, ranges, trimsection, medsection, outim,
+ index, c1, c2, l1, l2, isnull, median)
+
+pointer ir # pointer to the ir structure
+pointer imlist # pointer to the list of input images
+int ranges[ARB] # list of missing subrasters
+char trimsection[ARB] # input image section for output
+char medsection[ARB] # input image section for median computation
+pointer outim # pointer to the output image
+int index[ARB] # index array
+int c1[ARB] # array of beginning column limits
+int c2[ARB] # array of ending column limits
+int l1[ARB] # array of beginning line limits
+int l2[ARB] # array of ending line limits
+int isnull[ARB] # output input image order number
+real median[ARB] # output median of input image
+
+int i, j, k, nimrows, nimcols, imcount, next_null
+pointer sp, imname, im, buf
+int get_next_number(), imtgetim()
+pointer immap(), imgs2r()
+real amedr()
+
+begin
+ nimcols = IM_LEN(outim,1)
+ nimrows = IM_LEN(outim,2)
+
+ call smark (sp)
+ call salloc (imname, SZ_FNAME, TY_CHAR)
+
+ imcount = 1
+ next_null = 0
+ if (get_next_number (ranges, next_null) == EOF)
+ next_null = IR_NXSUB(ir) * IR_NYSUB(ir) + 1
+
+ # Loop over the input images.
+ do i = 1, IR_NXSUB(ir) * IR_NYSUB(ir) {
+
+ # Set the indices array.
+ call ir_indices (i, j, k, IR_NXSUB(ir), IR_NYSUB(ir),
+ IR_CORNER(ir), IR_RASTER(ir), IR_ORDER(ir))
+ index[i] = i
+ c1[i] = max (1, min (1 + (j - 1) * (IR_NCOLS(ir) -
+ IR_NXOVERLAP(ir)), nimcols))
+ c2[i] = min (nimcols, max (1, c1[i] + IR_NCOLS(ir) - 1))
+ l1[i] = max (1, min (1 + (k - 1) * (IR_NROWS(ir) -
+ IR_NYOVERLAP(ir)), nimrows))
+ l2[i] = min (nimrows, max (1, l1[i] + IR_NROWS(ir) - 1))
+
+ # Set the index of each image in the image template
+ # and compute the median of the subraster.
+ if (i < next_null) {
+ isnull[i] = imcount
+ if (medsection[1] != EOS) {
+ if (imtgetim (imlist, Memc[imname], SZ_FNAME) == EOF)
+ call error (0, "Error reading input image list.")
+ call strcat (medsection, Memc[imname], SZ_FNAME)
+ im = immap (Memc[imname], READ_ONLY, TY_CHAR)
+ buf = imgs2r (im, 1, int (IM_LEN(im,1)), 1, int (IM_LEN(im,
+ 2)))
+ median[i] = amedr (Memr[buf], int (IM_LEN(im,1)) *
+ int (IM_LEN(im,2)))
+ call imunmap (im)
+ } else
+ median[i] = INDEFR
+ imcount = imcount + 1
+ } else {
+ isnull[i] = 0
+ if (medsection[1] == EOS)
+ median[i] = INDEFR
+ else
+ median[i] = IR_OVAL(ir)
+ if (get_next_number (ranges, next_null) == EOF)
+ next_null = IR_NXSUB(ir) * IR_NYSUB(ir) + 1
+ }
+
+ }
+
+ call imtrew (imlist)
+ call sfree (sp)
+end
+
+
+# IR_SETIM -- Procedure to set up the output image characteristics.
+
+pointer procedure ir_setim (ir, list, trimsection, outimage, nimcols, nimrows,
+ opixtype)
+
+pointer ir # pointer to the ir structure
+pointer list # pointer to list of input images
+char trimsection[ARB]# input image section
+char outimage[ARB] # name of the output image
+int nimcols # number of output image columns
+int nimrows # number of output image rows
+int opixtype # output image pixel type
+
+int ijunk, nc, nr
+pointer sp, imname, im, outim
+pointer imtgetim(), immap()
+
+begin
+ call smark (sp)
+ call salloc (imname, SZ_FNAME, TY_CHAR)
+
+ # Get the size of the first subraster.
+ if (imtgetim (list, Memc[imname], SZ_FNAME) != EOF) {
+ call strcat (trimsection, Memc[imname], SZ_FNAME)
+ im = immap (Memc[imname], READ_ONLY, 0)
+ IR_NCOLS(ir) = IM_LEN(im,1)
+ IR_NROWS(ir) = IM_LEN(im,2)
+ call imunmap (im)
+ call imtrew (list)
+ } else
+ call error (0, "Error reading first input image.\n")
+
+ # Compute the size of the output image.
+ ijunk = IR_NXSUB(ir) * IR_NCOLS(ir) - (IR_NXSUB(ir) - 1) *
+ IR_NXOVERLAP(ir)
+ if (IS_INDEFI(nimcols))
+ nc = ijunk
+ else
+ nc = max (nimcols, ijunk)
+ ijunk = IR_NYSUB(ir) * IR_NROWS(ir) - (IR_NYSUB(ir) - 1) *
+ IR_NYOVERLAP(ir)
+ if (IS_INDEFI(ijunk))
+ nr = ijunk
+ else
+ nr = max (nimrows, ijunk)
+
+ # Set the output pixel type.
+ if (opixtype == ERR)
+ opixtype = TY_REAL
+
+ # Open output image and set the parameters.
+ outim = immap (outimage, NEW_IMAGE, 0)
+ IM_NDIM(outim) = 2
+ IM_LEN(outim,1) = nc
+ IM_LEN(outim,2) = nr
+ IM_PIXTYPE(outim) = opixtype
+
+ call sfree (sp)
+
+ return (outim)
+end
+
+
+# IR_MKMOSAIC -- Procedure to make the mosaiced image.
+
+procedure ir_mkmosaic (imlist, trimsection, outim, index, c1, c2, l1, l2,
+ isnull, median, nxsub, nysub, oval, subtract)
+
+pointer imlist # pointer to input image list
+char trimsection[ARB]# input image section
+pointer outim # pointer to the output image
+int index[ARB] # index array for sorting the images
+int c1[ARB] # array of column beginnings
+int c2[ARB] # array of column endings
+int l1[ARB] # array of line beginnings
+int l2[ARB] # array of line endings
+int isnull[ARB] # index of input image in the template
+real median[ARB] # array of input image median values
+int nxsub # number of subrasters per output image column
+int nysub # number of subrasters per output image row
+real oval # pixel value of undefined output image regions
+int subtract # subtract the median off each subraster
+
+int i, j, noutcols, noutlines, olineptr, ll1, ll2
+pointer sp, inimage, imptrs, buf
+pointer imtrgetim(), immap(), impl2r()
+
+begin
+ # Allocate temporary space.
+ call smark (sp)
+ call salloc (imptrs, nxsub, TY_POINTER)
+ call salloc (inimage, SZ_FNAME, TY_CHAR)
+
+ # Sort the subrasters on the yindex.
+ call ir_qsorti (l1, index, index, nxsub * nysub)
+
+ noutcols = IM_LEN(outim,1)
+ noutlines = IM_LEN(outim,2)
+
+ # Loop over the input images.
+ olineptr = 1
+ do i = 1, nxsub * nysub, nxsub {
+
+ # Compute the line and column limits.
+ ll1 = l1[index[i]]
+ ll2 = l2[index[i]]
+
+ # Open the nxsub input images.
+ do j = i, i + nxsub - 1 {
+ if (isnull[index[j]] <= 0) {
+ Memc[inimage] = EOS
+ Memi[imptrs+j-i] = NULL
+ } else {
+ if (imtrgetim (imlist, isnull[index[j]], Memc[inimage],
+ SZ_FNAME) == EOF)
+ Memi[imptrs+j-i] = NULL
+ else {
+ call strcat (trimsection, Memc[inimage], SZ_FNAME)
+ Memi[imptrs+j-i] = immap (Memc[inimage], READ_ONLY, 0)
+ }
+ }
+ }
+
+ # Write out the undefined lines.
+ while (olineptr < ll1) {
+ buf = impl2r (outim, olineptr)
+ call amovkr (oval, Memr[buf], noutcols)
+ olineptr = olineptr + 1
+ }
+
+ # Write the output lines.
+ call ir_mklines (Memi[imptrs], outim, index, c1, c2, ll1, ll2,
+ median, i, nxsub, oval, subtract)
+ olineptr = ll2 + 1
+
+ # Close up the images.
+ # Open the nxsub input images.
+ do j = i, i + nxsub - 1 {
+ if (Memi[imptrs+j-i] != NULL)
+ call imunmap (Memi[imptrs+j-i])
+ }
+
+ }
+
+ # Write out the remaining undefined lines.
+ while (olineptr < noutlines) {
+ buf = impl2r (outim, olineptr)
+ call amovkr (oval, Memr[buf], noutcols)
+ olineptr = olineptr + 1
+ }
+
+ call sfree (sp)
+end
+
+
+# IR_MKLINES -- Construct and output image lines.
+
+procedure ir_mklines (imptrs, outim, index, c1, c2, l1, l2, meds, init, nsub,
+ oval, subtract)
+
+pointer imptrs[ARB] # array of input image pointers
+pointer outim # output imnage pointer
+int index[ARB] # array of indices
+int c1[ARB] # array of beginning columns
+int c2[ARB] # array of ending columns
+int l1 # beginning line
+int l2 # ending line
+real meds[ARB] # array of median values
+int init # first index
+int nsub # number of subrasters
+real oval # output value
+int subtract # subtract the median value
+
+int i, j, jj, noutcols
+pointer obuf, ibuf
+pointer impl2r(), imgl2r()
+
+begin
+ noutcols = IM_LEN(outim, 1)
+ do i = l1, l2 {
+ obuf = impl2r (outim, i)
+ call amovkr (oval, Memr[obuf], noutcols)
+ do j = 1, nsub {
+ jj = index[j+init-1]
+ if (imptrs[j] != NULL) {
+ ibuf = imgl2r (imptrs[j], i - l1 + 1)
+ if (subtract == YES)
+ call asubkr (Memr[ibuf], meds[jj], Memr[obuf+c1[jj]-1],
+ c2[jj] - c1[jj] + 1)
+ else
+ call amovr (Memr[ibuf], Memr[obuf+c1[jj]-1], c2[jj] -
+ c1[jj] + 1)
+ }
+ }
+ }
+end
+
+
+# IR_DTWINPUT -- Procedure to write the output database file.
+
+procedure ir_dtwinput (imlist, trimsection, outimage, dt, index, c1, c2, l1,
+ l2, isnull, median, nsub, subtract, verbose)
+
+int imlist # input image list
+char trimsection[ARB]# trim section of input image
+char outimage[ARB] # output image
+pointer dt # pointer to the database file
+int index[ARB] # array of sorted indices (not used at present)
+int c1[ARB] # array of beginning column limits
+int c2[ARB] # array of ending column limits
+int l1[ARB] # array of beginning line limits
+int l2[ARB] # array of ending line limits
+int isnull[ARB] # image name index
+real median[ARB] # array of medians
+int nsub # number of subrasters
+int subtract # subtract the median from the subraster
+int verbose # print verbose messages
+
+int i
+pointer sp, imname
+int imtrgetim()
+
+begin
+ call smark (sp)
+ call salloc (imname, SZ_FNAME, TY_CHAR)
+
+ # Write out the number of subrasters.
+ call dtput (dt, "\tnsubrasters\t%d\n")
+ call pargi (nsub)
+
+ do i = 1, nsub {
+
+ if (isnull[i] <= 0)
+ call strcpy ("nullimage", Memc[imname], SZ_FNAME)
+ else if (imtrgetim (imlist, isnull[i], Memc[imname],
+ SZ_FNAME) != EOF)
+ call strcat (trimsection, Memc[imname], SZ_FNAME)
+ else
+ Memc[imname] = EOS
+
+ call dtput (dt,"\t%s %s[%d:%d,%d:%d] %g %g\n")
+ call pargstr (Memc[imname])
+ call pargstr (outimage)
+ call pargi (c1[i])
+ call pargi (c2[i])
+ call pargi (l1[i])
+ call pargi (l2[i])
+ call pargr (median[i])
+ if (subtract == YES)
+ call pargr (-median[i])
+ else
+ call pargr (0.0)
+
+ if (verbose == YES) {
+ call printf ("imcopy %s %s[%d:%d,%d:%d] %g %g\n")
+ call pargstr (Memc[imname])
+ call pargstr (outimage)
+ call pargi (c1[i])
+ call pargi (c2[i])
+ call pargi (l1[i])
+ call pargi (l2[i])
+ call pargr (median[i])
+ if (subtract == YES)
+ call pargr (-median[i])
+ else
+ call pargr (0.0)
+ }
+ }
+
+ call sfree (sp)
+end
generated by cgit (git 2.34.1) at 2025-09-20 16:21:41 -0400