diff options
Diffstat (limited to 'noao/nproto/ir')
| -rw-r--r-- | noao/nproto/ir/iralign.h | 55 | ||||
| -rw-r--r-- | noao/nproto/ir/iralign.x | 376 | ||||
| -rw-r--r-- | noao/nproto/ir/irdbio.x | 117 | ||||
| -rw-r--r-- | noao/nproto/ir/iriinit.x | 28 | ||||
| -rw-r--r-- | noao/nproto/ir/irimisec.x | 105 | ||||
| -rw-r--r-- | noao/nproto/ir/irimzero.x | 22 | ||||
| -rw-r--r-- | noao/nproto/ir/irindices.x | 139 | ||||
| -rw-r--r-- | noao/nproto/ir/irlinks.x | 496 | ||||
| -rw-r--r-- | noao/nproto/ir/irmatch1d.x | 122 | ||||
| -rw-r--r-- | noao/nproto/ir/irmatch2d.x | 276 | ||||
| -rw-r--r-- | noao/nproto/ir/irmedr.x | 35 | ||||
| -rw-r--r-- | noao/nproto/ir/iroverlap.x | 40 | ||||
| -rw-r--r-- | noao/nproto/ir/irqsort.x | 215 | ||||
| -rw-r--r-- | noao/nproto/ir/irtools.x | 147 | ||||
| -rw-r--r-- | noao/nproto/ir/mkpkg | 24 | ||||
| -rw-r--r-- | noao/nproto/ir/t_iralign.x | 134 | ||||
| -rw-r--r-- | noao/nproto/ir/t_irmatch1d.x | 159 | ||||
| -rw-r--r-- | noao/nproto/ir/t_irmatch2d.x | 159 | ||||
| -rw-r--r-- | noao/nproto/ir/t_irmosaic.x | 498 |
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 |