From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- pkg/xtools/fixpix/xtpmmap.x | 693 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 693 insertions(+) create mode 100644 pkg/xtools/fixpix/xtpmmap.x (limited to 'pkg/xtools/fixpix/xtpmmap.x') diff --git a/pkg/xtools/fixpix/xtpmmap.x b/pkg/xtools/fixpix/xtpmmap.x new file mode 100644 index 00000000..54bbf954 --- /dev/null +++ b/pkg/xtools/fixpix/xtpmmap.x @@ -0,0 +1,693 @@ +include +include +include +include +include +include +include +include + + +# XT_PMMAP -- Open a pixel mask READ_ONLY. +# +# This routine maps multiple types of mask files and designations. +# It matches the mask coordinates to the reference image based on the +# physical coordinate system so the mask may be of a different size. +# The mask name is returned so that the task has the name pointed to by "BPM". +# A null filename is allowed and returns NULL. + +pointer procedure xt_pmmap (pmname, refim, mname, sz_mname) + +char pmname[ARB] #I Pixel mask name +pointer refim #I Reference image pointer +char mname[ARB] #O Expanded mask name +int sz_mname #O Size of expanded mask name + +int i, flag, nowhite() +pointer sp, fname, im, ref, xt_pmmap1() +bool streq() +errchk xt_pmmap1 + +begin + call smark (sp) + call salloc (fname, SZ_FNAME, TY_CHAR) + + im = NULL + i = nowhite (pmname, Memc[fname], SZ_FNAME) + if (Memc[fname] == '!') { + iferr (call imgstr (refim, Memc[fname+1], Memc[fname], SZ_FNAME)) + Memc[fname] = EOS + } else if (streq (Memc[fname], "BPM")) { + iferr (call imgstr (refim, "BPM", Memc[fname], SZ_FNAME)) + Memc[fname] = EOS + } else if (streq (Memc[fname], "^BPM")) { + flag = INVERT_MASK + iferr (call imgstr (refim, "BPM", Memc[fname+1], SZ_FNAME)) + Memc[fname] = EOS + } + + if (Memc[fname] == '^') { + flag = INVERT_MASK + call strcpy (Memc[fname+1], Memc[fname], SZ_FNAME) + } else + flag = NO + + if (streq (Memc[fname], "EMPTY")) + ref = refim + else + ref = NULL + + if (Memc[fname] != EOS) + im = xt_pmmap1 (Memc[fname], ref, refim, flag, YES) + call strcpy (Memc[fname], mname, sz_mname) + + call sfree (sp) + return (im) +end + + +# XT_MAPPM -- Open a pixel mask READ_ONLY with/without matching. +# +# This routine maps multiple types of mask files and designations. +# It may match the mask coordinates to the reference image based on the +# physical coordinate system. In either case the mask is matched to be +# the same size. The mask name is returned so that the task has the +# name pointed to by "BPM". A null filename is allowed and returns NULL. + +pointer procedure xt_mappm (pmname, refim, match, mname, sz_mname) + +char pmname[ARB] #I Pixel mask name +pointer refim #I Reference image pointer +int match #I Match by physical coordinates? +char mname[ARB] #O Expanded mask name +int sz_mname #O Size of expanded mask name + +int i, flag, nowhite() +pointer sp, fname, im, ref, xt_pmmap1() +bool streq() +errchk xt_pmmap1 + +begin + call smark (sp) + call salloc (fname, SZ_FNAME, TY_CHAR) + + im = NULL + i = nowhite (pmname, Memc[fname], SZ_FNAME) + if (Memc[fname] == '!') { + iferr (call imgstr (refim, Memc[fname+1], Memc[fname], SZ_FNAME)) + Memc[fname] = EOS + } else if (streq (Memc[fname], "BPM")) { + iferr (call imgstr (refim, "BPM", Memc[fname], SZ_FNAME)) + Memc[fname] = EOS + } else if (streq (Memc[fname], "^BPM")) { + flag = INVERT_MASK + iferr (call imgstr (refim, "BPM", Memc[fname+1], SZ_FNAME)) + Memc[fname] = EOS + } + + if (Memc[fname] == '^') { + flag = INVERT_MASK + call strcpy (Memc[fname+1], Memc[fname], SZ_FNAME) + } else + flag = NO + + if (streq (Memc[fname], "EMPTY")) + ref = refim + else + ref = NULL + + if (Memc[fname] != EOS) + im = xt_pmmap1 (Memc[fname], ref, refim, flag, match) + call strcpy (Memc[fname], mname, sz_mname) + + call sfree (sp) + return (im) +end + + +# XT_PMUNMAP -- Unmap a mask image. +# Note that the imio pointer may be purely an internal pointer opened +# with im_pmmapo so we need to free the pl pointer explicitly. + +procedure xt_pmunmap (im) + +pointer im #I IMIO pointer for mask + +pointer pm +int imstati() + +begin + pm = imstati (im, IM_PMDES) + call pm_close (pm) + call imseti (im, IM_PMDES, NULL) + call imunmap (im) +end + + +# XT_PMMAP1 -- Open a pixel mask READ_ONLY. The input mask may be +# a pixel list image, a non-pixel list image, or a text file. +# Return error if the pixel mask cannot be opened. For pixel masks +# or image masks match the WCS. + +pointer procedure xt_pmmap1 (pmname, ref, refim, flag, match) + +char pmname[ARB] #I Pixel mask name +pointer ref #I Reference image for pixel mask +pointer refim #I Reference image for image or text +int flag #I Mask flag +int match #I Match by physical coordinates? + +int imstati(), errcode() +pointer im, pm +pointer im_pmmap(), xt_pmimmap(), xt_pmtext(), xt_pmsection() +bool streq() +errchk xt_match + +begin + im = NULL + + if (streq (pmname, "STDIN")) + im = xt_pmtext (pmname, refim, flag) + + else if (pmname[1] == '[') + im = xt_pmsection (pmname, refim, flag) + + else { + ifnoerr (im = im_pmmap (pmname, READ_ONLY, ref)) { + call xt_match (im, refim, match) + if (flag == INVERT_MASK) { + pm = imstati (im, IM_PMDES) + call xt_pminvert (pm) + call imseti (im, IM_PMDES, pm) + } + } else { + switch (errcode()) { + case SYS_IKIOPEN, SYS_FOPNNEXFIL, SYS_PLBADSAVEF, SYS_FOPEN: + ifnoerr (im = xt_pmimmap (pmname, refim, flag)) + call xt_match (im, refim, match) + else { + switch (errcode()) { + case SYS_IKIOPEN: + im = xt_pmtext (pmname, refim, flag) + default: + call erract (EA_ERROR) + } + } + default: + call erract (EA_ERROR) + } + } + } + + return (im) +end + + +# XT_PMIMMAP -- Open a pixel mask from a non-pixel list image. +# Return error if the image cannot be opened. + +pointer procedure xt_pmimmap (pmname, refim, flag) + +char pmname[ARB] #I Image name +pointer refim #I Reference image pointer +int flag #I Mask flag + +int i, ndim, npix, rop, val +pointer sp, v1, v2, im_in, im_out, pm, mw, data + +int imstati(), imgnli() +pointer immap(), pm_newmask(), im_pmmapo(), imgl1i(), mw_openim() +errchk immap, mw_openim, im_pmmapo + +begin + call smark (sp) + call salloc (v1, IM_MAXDIM, TY_LONG) + call salloc (v2, IM_MAXDIM, TY_LONG) + + call amovkl (long(1), Meml[v1], IM_MAXDIM) + call amovkl (long(1), Meml[v2], IM_MAXDIM) + + im_in = immap (pmname, READ_ONLY, 0) + pm = imstati (im_in, IM_PMDES) + if (pm != NULL) + return (im_in) + pm = pm_newmask (im_in, 16) + + ndim = IM_NDIM(im_in) + npix = IM_LEN(im_in,1) + + if (flag == INVERT_MASK) + rop = PIX_NOT(PIX_SRC) + else + rop = PIX_SRC + + while (imgnli (im_in, data, Meml[v1]) != EOF) { + if (flag == INVERT_MASK) { + do i = 0, npix-1 { + val = Memi[data+i] + if (val <= 0) + Memi[data+i] = 1 + else + Memi[data+i] = 0 + } + } else { + do i = 0, npix-1 { + val = Memi[data+i] + if (val < 0) + Memi[data+i] = 0 + } + } + call pmplpi (pm, Meml[v2], Memi[data], 0, npix, rop) + call amovl (Meml[v1], Meml[v2], ndim) + } + + im_out = im_pmmapo (pm, im_in) + data = imgl1i (im_out) # Force I/O to set header + mw = mw_openim (im_in) # Set WCS + call mw_saveim (mw, im_out) + call mw_close (mw) + + #call imunmap (im_in) + call xt_pmunmap (im_in) + call sfree (sp) + return (im_out) +end + + +# XT_PMTEXT -- Create a pixel mask from a text file of rectangles. +# Return error if the file cannot be opened. +# This routine only applies to the first 2D plane. + +pointer procedure xt_pmtext (pmname, refim, flag) + +char pmname[ARB] #I Image name +pointer refim #I Reference image pointer +int flag #I Mask flag + +int fd, nc, nl, c1, c2, l1, l2, nc1, nl1, rop +pointer pm, im, mw, dummy + +int open(), fscan(), nscan() +pointer pm_newmask(), im_pmmapo(), imgl1i(), mw_openim() +errchk open, im_pmmapo + +begin + fd = open (pmname, READ_ONLY, TEXT_FILE) + pm = pm_newmask (refim, 16) + + nc = IM_LEN(refim,1) + nl = IM_LEN(refim,2) + + if (flag == INVERT_MASK) + call pl_box (pm, 1, 1, nc, nl, PIX_SET+PIX_VALUE(1)) + + while (fscan (fd) != EOF) { + call gargi (c1) + call gargi (c2) + call gargi (l1) + call gargi (l2) + if (nscan() != 4) { + if (nscan() == 2) { + l1 = c2 + c2 = c1 + l2 = l1 + } else + next + } + + c1 = max (1, c1) + c2 = min (nc, c2) + l1 = max (1, l1) + l2 = min (nl, l2) + nc1 = c2 - c1 + 1 + nl1 = l2 - l1 + 1 + if (nc1 < 1 || nl1 < 1) + next + + # Select mask value based on shape of rectangle. + if (flag == INVERT_MASK) + rop = PIX_CLR + else if (nc1 <= nl1) + rop = PIX_SET+PIX_VALUE(2) + else + rop = PIX_SET+PIX_VALUE(3) + + # Set mask rectangle. + call pm_box (pm, c1, l1, c2, l2, rop) + } + + call close (fd) + im = im_pmmapo (pm, refim) + dummy = imgl1i (im) # Force I/O to set header + mw = mw_openim (refim) # Set WCS + call mw_saveim (mw, im) + call mw_close (mw) + + return (im) +end + + +# XT_PMSECTION -- Create a pixel mask from an image section. +# This only applies the mask to the first plane of the image. + +pointer procedure xt_pmsection (section, refim, flag) + +char section[ARB] #I Image section +pointer refim #I Reference image pointer +int flag #I Mask flag + +int i, j, ip, temp, a[2], b[2], c[2], rop, ctoi() +pointer pm, im, mw, dummy, pm_newmask(), im_pmmapo(), imgl1i(), mw_openim() +errchk im_pmmapo +define error_ 99 + +begin + # This is currently only for 1D and 2D images. + if (IM_NDIM(refim) > 2) + call error (1, "Image sections only allowed for 1D and 2D images") + + # Decode the section string. + call amovki (1, a, 2) + call amovki (1, b, 2) + call amovki (1, c, 2) + do i = 1, IM_NDIM(refim) + b[i] = IM_LEN(refim,i) + + ip = 1 + while (IS_WHITE(section[ip])) + ip = ip + 1 + if (section[ip] == '[') { + ip = ip + 1 + + do i = 1, IM_NDIM(refim) { + while (IS_WHITE(section[ip])) + ip = ip + 1 + + # Get a:b:c. Allow notation such as "-*:c" + # (or even "-:c") where the step is obviously negative. + + if (ctoi (section, ip, temp) > 0) { # a + a[i] = temp + if (section[ip] == ':') { + ip = ip + 1 + if (ctoi (section, ip, b[i]) == 0) # a:b + goto error_ + } else + b[i] = a[i] + } else if (section[ip] == '-') { # -* + temp = a[i] + a[i] = b[i] + b[i] = temp + ip = ip + 1 + if (section[ip] == '*') + ip = ip + 1 + } else if (section[ip] == '*') # * + ip = ip + 1 + if (section[ip] == ':') { # ..:step + ip = ip + 1 + if (ctoi (section, ip, c[i]) == 0) + goto error_ + else if (c[i] == 0) + goto error_ + } + if (a[i] > b[i] && c[i] > 0) + c[i] = -c[i] + + while (IS_WHITE(section[ip])) + ip = ip + 1 + if (i < IM_NDIM(refim)) { + if (section[ip] != ',') + goto error_ + } else { + if (section[ip] != ']') + goto error_ + } + ip = ip + 1 + } + } + + # In this case make the values be increasing only. + do i = 1, IM_NDIM(refim) + if (c[i] < 0) { + temp = a[i] + a[i] = b[i] + b[i] = temp + c[i] = -c[i] + } + + # Make the mask. + pm = pm_newmask (refim, 16) + + if (flag == INVERT_MASK) { + rop = PIX_SET+PIX_VALUE(1) + call pm_box (pm, 1, 1, IM_LEN(refim,1), IM_LEN(refim,2), rop) + rop = PIX_CLR + } else + rop = PIX_SET+PIX_VALUE(1) + + if (c[1] == 1 && c[2] == 1) + call pm_box (pm, a[1], a[2], b[1], b[2], rop) + + else if (c[1] == 1) + for (i=a[2]; i<=b[2]; i=i+c[2]) + call pm_box (pm, a[1], i, b[1], i, rop) + + else + for (i=a[2]; i<=b[2]; i=i+c[2]) + for (j=a[1]; j<=b[1]; j=j+c[1]) + call pm_point (pm, j, i, rop) + + im = im_pmmapo (pm, refim) + dummy = imgl1i (im) # Force I/O to set header + mw = mw_openim (refim) # Set WCS + call mw_saveim (mw, im) + call mw_close (mw) + + return (im) + +error_ + call error (1, "Error in image section specification") +end + + +# XT_PMINVERT -- Invert a pixel mask by changing 0 to 1 and non-zero to zero. + +procedure xt_pminvert (pm) + +pointer pm #I Pixel mask to be inverted + +int i, naxes, axlen[IM_MAXDIM], depth, npix, val +pointer sp, v, buf, one +bool pm_linenotempty() + +begin + call pm_gsize (pm, naxes, axlen, depth) + + call smark (sp) + call salloc (v, IM_MAXDIM, TY_LONG) + call salloc (buf, axlen[1], TY_INT) + call salloc (one, 6, TY_INT) + + npix = axlen[1] + RLI_LEN(one) = 2 + RLI_AXLEN(one) = npix + Memi[one+3] = 1 + Memi[one+4] = npix + Memi[one+5] = 1 + + call amovkl (long(1), Meml[v], IM_MAXDIM) + repeat { + if (pm_linenotempty (pm, Meml[v])) { + call pmglpi (pm, Meml[v], Memi[buf], 0, npix, 0) + do i = 0, npix-1 { + val = Memi[buf+i] + if (val == 0) + Memi[buf+i] = 1 + else + Memi[buf+i] = 0 + } + call pmplpi (pm, Meml[v], Memi[buf], 0, npix, PIX_SRC) + } else + call pmplri (pm, Meml[v], Memi[one], 0, npix, PIX_SRC) + + do i = 2, naxes { + Meml[v+i-1] = Meml[v+i-1] + 1 + if (Meml[v+i-1] <= axlen[i]) + break + else if (i < naxes) + Meml[v+i-1] = 1 + } + } until (Meml[v+naxes-1] > axlen[naxes]) + + call sfree (sp) +end + + +# XT_MATCH -- Set the pixel mask to match the reference image. +# This matches sizes and physical coordinates and allows the +# original mask to be smaller or larger than the reference image. +# Subsequent use of the pixel mask can then work in the logical +# coordinates of the reference image. The mask values are the maximum +# of the mask values which overlap each reference image pixel. +# A null input returns a null output. + +procedure xt_match (im, refim, match) + +pointer im #U Pixel mask image pointer +pointer refim #I Reference image pointer +int match #I Match by physical coordinates? + +int i, j, k, l, i1, i2, j1, j2, nc, nl, ncpm, nlpm, nx, val +double x1, x2, y1, y2, lt[6], lt1[6], lt2[6] +long vold[IM_MAXDIM], vnew[IM_MAXDIM] +pointer pm, pmnew, imnew, mw, ctx, cty, bufref, bufpm + +int imstati() +pointer pm_open(), mw_openim(), im_pmmapo(), imgl1i(), mw_sctran() +bool pm_empty(), pm_linenotempty() +errchk pm_open, mw_openim, im_pmmapo + +begin + if (im == NULL) + return + + # Set sizes. + nc = IM_LEN(refim,1) + nl = IM_LEN(refim,2) + ncpm = IM_LEN(im,1) + nlpm = IM_LEN(im,2) + + # If the mask is empty and the sizes are the same then it does not + # matter if the two are actually matched in physical coordinates. + pm = imstati (im, IM_PMDES) + if (pm_empty(pm) && nc == ncpm && nl == nlpm) + return + + # Compute transformation between reference (logical) coordinates + # and mask (physical) coordinates if desired. + + mw = mw_openim (im) + call mw_gltermd (mw, lt, lt[5], 2) + call mw_close (mw) + + if (match == YES) { + mw = mw_openim (refim) + call mw_gltermd (mw, lt2, lt2[5], 2) + call mw_close (mw) + } else + call amovd (lt, lt2, 6) + + # Combine lterms. + call mw_invertd (lt, lt1, 2) + call mw_mmuld (lt1, lt2, lt, 2) + call mw_vmuld (lt, lt[5], lt[5], 2) + lt[5] = lt2[5] - lt[5] + lt[6] = lt2[6] - lt[6] + do i = 1, 6 + lt[i] = nint (1D6 * (lt[i]-int(lt[i]))) / 1D6 + int(lt[i]) + + # Check for a rotation. For now don't allow any rotation. + if (lt[2] != 0. || lt[3] != 0.) + call error (1, "Image and mask have a relative rotation") + + # Check for an exact match. + if (lt[1] == 1D0 && lt[4] == 1D0 && lt[5] == 0D0 && lt[6] == 0D0 && + nc == ncpm && nl == nlpm) + return + + # Set reference to mask coordinates. + mw = mw_openim (im) + call mw_sltermd (mw, lt, lt[5], 2) + ctx = mw_sctran (mw, "logical", "physical", 1) + cty = mw_sctran (mw, "logical", "physical", 2) + + # Create a new pixel mask of the required size and offset. + # Do dummy image I/O to set the header. + pmnew = pm_open (NULL) + call pm_ssize (pmnew, 2, IM_LEN(refim,1), 27) + imnew = im_pmmapo (pmnew, NULL) + bufref = imgl1i (imnew) + + # Compute region of mask overlapping the reference image. + call mw_ctrand (ctx, 1-0.5D0, x1, 1) + call mw_ctrand (ctx, nc+0.5D0, x2, 1) + i1 = max (1, nint(min(x1,x2)+1D-5)) + i2 = min (ncpm, nint(max(x1,x2)-1D-5)) + call mw_ctrand (cty, 1-0.5D0, y1, 1) + call mw_ctrand (cty, nl+0.5D0, y2, 1) + j1 = max (1, nint(min(y1,y2)+1D-5)) + j2 = min (nlpm, nint(max(y1,y2)-1D-5)) + + # Set the new mask values to the maximum of all mask values falling + # within each reference pixel in the overlap region. + if (i1 <= i2 && j1 <= j2) { + nx = i2 - i1 + 1 + vold[1] = i1 + vnew[1] = 1 + + # If the scales are the same then it is just a problem of + # padding. In this case use range lists for speed. + if (lt[1] == 1D0 && lt[4] == 1D0) { + call malloc (bufpm, 3+3*nc, TY_INT) + k = nint (lt[5]) + l = nint (lt[6]) + do j = max(1-l,j1), min(nl-l,j2) { + vold[2] = j + call pmglri (pm, vold, Memi[bufpm], 0, nc, PIX_SRC) + if (k != 0) { + bufref = bufpm + do i = 2, Memi[bufpm] { + bufref = bufref + 3 + Memi[bufref] = Memi[bufref] + k + } + } + vnew[2] = j + l + call pmplri (pmnew, vnew, Memi[bufpm], 0, nc, PIX_SRC) + } + bufref = NULL + + # Do all the geometry and pixel size matching. This can + # be slow. + } else { + call malloc (bufpm, nx, TY_INT) + call malloc (bufref, nc, TY_INT) + do j = 1, nl { + call mw_ctrand (cty, j-0.5D0, y1, 1) + call mw_ctrand (cty, j+0.5D0, y2, 1) + j1 = max (1, nint(min(y1,y2)+1D-5)) + j2 = min (nlpm, nint(max(y1,y2)-1D-5)) + if (j2 < j1) + next + + vnew[2] = j + call aclri (Memi[bufref], nc) + do l = j1, j2 { + vold[2] = l + if (!pm_linenotempty (pm, vold)) + next + call pmglpi (pm, vold, Memi[bufpm], 0, nx, 0) + do i = 1, nc { + call mw_ctrand (ctx, i-0.5D0, x1, 1) + call mw_ctrand (ctx, i+0.5D0, x2, 1) + i1 = max (1, nint(min(x1,x2)+1D-5)) + i2 = min (ncpm, nint(max(x1,x2)-1D-5)) + if (i2 < i1) + next + val = Memi[bufref+i-1] + do k = i1-vold[1], i2-vold[1] + val = max (val, Memi[bufpm+k]) + Memi[bufref+i-1] = val + } + } + call pmplpi (pmnew, vnew, Memi[bufref], 0, nc, PIX_SRC) + } + } + call mfree (bufref, TY_INT) + call mfree (bufpm, TY_INT) + } + + call mw_close (mw) + call xt_pmunmap (im) + im = imnew + call imseti (im, IM_PMDES, pmnew) +end -- cgit