diff options
Diffstat (limited to 'noao/nproto/ace/skyblock.x')
| -rw-r--r-- | noao/nproto/ace/skyblock.x | 1039 |
1 files changed, 1039 insertions, 0 deletions
diff --git a/noao/nproto/ace/skyblock.x b/noao/nproto/ace/skyblock.x new file mode 100644 index 00000000..5e3eb5f9 --- /dev/null +++ b/noao/nproto/ace/skyblock.x @@ -0,0 +1,1039 @@ +include <error.h> +include <ctype.h> +include <imhdr.h> +include <imset.h> +include <mach.h> +include "skyblock.h" + + +# SKY_BLOCK - Determine sky and sky sigma in blocks. +# +# This is layered on MAPIO and CONVOLVE. + +procedure sky_block (skb, dosky, dosig, in, bpm, expmap, skyname, signame, + skymap, sigmap, logfd) + +pointer skb #U Sky block structure +bool dosky #I Compute sky +bool dosig #I Compute sigma +pointer in #I Input image pointer +pointer bpm #I Input mask +pointer expmap #I Exposure map +char skyname[ARB] #I Sky map name (if none then no output) +char signame[ARB] #I Sigma map name (if none then no output) +pointer skymap #U Sky map +pointer sigmap #U Sigma map +int logfd #I Verbose? + +int l, blkstep, nc, nl +real cnvwt +pointer sp, cnv, cnvdata, bp +pointer im[2], indata, skydata, sigdata, expdata +errchk skb_pars, skb_iminit, convolve, skb_accum, skb_update + +begin + if (!(dosky||dosig)) + return + + call smark (sp) + + # Log operation. + if (logfd != NULL) { + if (dosky && dosig) + call fprintf (logfd, + " Determine sky and sigma by block statistics:\n") + else if (dosky) + call fprintf (logfd, " Determine sky by block statistics:\n") + else + call fprintf (logfd, + " Determine sigma by block statistics:\n") + } + + # Set parameters if not set in a previous call or set externally. + if (skb == NULL) + call skb_pars ("open", "", skb) + + # Set parameters for the image. + blkstep = SKB_BLKSTEP(skb) + call skb_iminit (skb, in, expmap, blkstep, logfd) + + # Set maximum number of image columns and lines to use. + nc = SKB_NCSBLK(skb) * SKB_NCSPIX(skb) + nl = SKB_NLSBLK(skb) * SKB_NLSPIX(skb) + + # Set up convolution. Note we can't use convolution with a blkstep. + cnv = SKB_CNV(skb) + if (Memc[cnv] != EOS) { + if (blkstep > 1) { + call salloc (cnv, 1, TY_CHAR) + Memc[cnv] = EOS + } else + call salloc (cnvdata, nc, TY_REAL) + } + + # Setup bad pixel mask. + if (bpm == NULL) { + call salloc (bp, nc, TY_INT) + call aclri (Memi[bp], nc) + } + + # Go through image creating low resolution sky blocks. + im[1] = in; im[2] = NULL + do l = 1, nl, blkstep { + call convolve (im, bpm, skymap, sigmap, expmap, 0, + 1., l, Memc[cnv], indata, bp, cnvdata, skydata, + sigdata, expdata, cnvwt, logfd) + call skb_accum (skb, l, blkstep, Memr[cnvdata], Memr[skydata], + Memr[sigdata], Memr[expdata], Memi[bp], nc, cnvwt) + } + + # Free convolution memory. + call convolve (im, bpm, skymap, sigmap, expmap, 0, + 1., 0, Memc[cnv], indata, bp, cnvdata, skydata, + sigdata, expdata, cnvwt, logfd) + + # Turn the sky blocks into sky maps. + call skb_update (skb, dosky, dosig, in, skyname, signame, + skymap, sigmap, logfd) + + # Free memory. + call skb_imfree (skb) + call sfree (sp) +end + + +# SKB_IMINIT -- Initialize parameters and allocate memory for an image. + +procedure skb_iminit (skb, im, expmap, blkstep, logfd) + +pointer skb #U Sky block structure +pointer im #I Image pointer +pointer expmap #I Exposure map pointer +int blkstep #U Line step for speed +int logfd #I Log file descriptor + +int nc, nl + +begin + # Number of pixels per subblock. + nc = IM_LEN(im,1) + nl = IM_LEN(im,2) + if (SKB_BLKSIZE(skb) < 0) { + if (nc < nl) { + SKB_NCSPIX(skb) = max (SKB_NMINPIX(skb), + nc / (SKB_NSUBBLKS(skb) * max(1,-SKB_BLKSIZE(skb)))) + SKB_NLSPIX(skb) = SKB_NCSPIX(skb) + } else { + SKB_NLSPIX(skb) = max (SKB_NMINPIX(skb), + nl / (SKB_NSUBBLKS(skb) * max(1,-SKB_BLKSIZE(skb)))) + SKB_NCSPIX(skb) = SKB_NLSPIX(skb) + } + } else { + SKB_NCSPIX(skb) = max (SKB_NMINPIX(skb), + min (nc, SKB_BLKSIZE(skb)) / SKB_NSUBBLKS(skb)) + SKB_NLSPIX(skb) = max (SKB_NMINPIX(skb), + min (nl, SKB_BLKSIZE(skb)) / SKB_NSUBBLKS(skb)) + } + + # Number of subblocks, blocks, and number of pixels per block. + SKB_NCSBLK(skb) = max (1, nc / SKB_NCSPIX(skb)) + SKB_NLSBLK(skb) = max (1, nl / SKB_NLSPIX(skb)) + SKB_NCBLK(skb) = (SKB_NCSBLK(skb)+SKB_NSUBBLKS(skb)-1)/SKB_NSUBBLKS(skb) + SKB_NLBLK(skb) = (SKB_NLSBLK(skb)+SKB_NSUBBLKS(skb)-1)/SKB_NSUBBLKS(skb) + SKB_NCPIX(skb) = SKB_NCSPIX(skb) * SKB_NSUBBLKS(skb) + SKB_NLPIX(skb) = SKB_NLSPIX(skb) * SKB_NSUBBLKS(skb) + + # Each subblock must have at least SKYMIN or FRAC sky pixels. + SKB_NSKYMIN(skb) = min (SKB_SKYMIN(skb), + nint (SKB_FRAC(skb) * SKB_NCSPIX(skb) * SKB_NLSPIX(skb))) + + # Histogram parameters. + SKB_NAV(skb) = nint (real(SKB_NBINS(skb)) / (min (SKB_NBINS(skb), + SKB_NCSPIX(skb) * SKB_NLSPIX(skb) / SKB_NMINBINS(skb)))) + SKB_NAV(skb) = SKB_NAV(skb) + mod (SKB_NAV(skb)+1, 2) + #SKB_NAV(skb) = 1 + + # Set line subsampling for speed. + if (blkstep > 1) { + blkstep = min (1 + SKB_NLSPIX(skb) / 30, blkstep) + SKB_NSKYMIN(skb) = SKB_NSKYMIN(skb) / blkstep + } + + # Allocate and initialize memory. + call calloc (SKB_BINS(skb), SKB_NBINS(skb)*(SKB_NCSBLK(skb)+1), TY_INT) + call calloc (SKB_NSKY(skb), SKB_NCSBLK(skb), TY_INT) + call calloc (SKB_SKYS(skb), SKB_NCSBLK(skb)*SKB_NLSBLK(skb), TY_REAL) + call calloc (SKB_SIGS(skb), SKB_NCSBLK(skb)*SKB_NLSBLK(skb), TY_REAL) + if (expmap == NULL) { + call malloc (SKB_EXP(skb), 1, TY_REAL) + Memr[SKB_EXP(skb)] = INDEFR + } else + call calloc (SKB_EXP(skb), SKB_NCSBLK(skb), TY_REAL) + + # Set pointers to first line of blocks. + SKB_SKY(skb) = SKB_SKYS(skb) + SKB_SIG(skb) = SKB_SIGS(skb) + + if (logfd != NULL) { + call fprintf (logfd, " Number of blocks: %d %d\n") + call pargi (SKB_NCBLK(skb)) + call pargi (SKB_NLBLK(skb)) + call fprintf (logfd, " Number of pixels per block: %d %d\n") + call pargi (SKB_NCPIX(skb)) + call pargi (SKB_NLPIX(skb)) + call fprintf (logfd, " Number of subblocks: %d %d\n") + call pargi (SKB_NCSBLK(skb)) + call pargi (SKB_NLSBLK(skb)) + call fprintf (logfd, " Number of pixels per subblock: %d %d\n") + call pargi (SKB_NCSPIX(skb)) + call pargi (SKB_NLSPIX(skb)) + if (blkstep > 1) { + call fprintf (logfd, " Line sampling step: %d\n") + call pargi (blkstep) + } + } +end + + +# SKB_IMFREE -- Free memory for an image. + +procedure skb_imfree (skb) + +pointer skb #I Sky block structure + +begin + call mfree (SKB_BINS(skb), TY_INT) + call mfree (SKB_NSKY(skb), TY_INT) + call mfree (SKB_SKYS(skb), TY_REAL) + call mfree (SKB_SIGS(skb), TY_REAL) + call mfree (SKB_EXP(skb), TY_REAL) +end + + +# SKB_ACCUM -- Accumulate sky pixels in block histograms. +# Evaluate histograms when a block is complete. + +procedure skb_accum (skb, line, blkstep, cnv, sky, sig, exp, bp, nc, cnvwt) + +pointer skb #I Sky block structure +int line #I Line +int blkstep #I Line step +real cnv[nc] #I Convolved image data +real sky[nc] #I Sky data +real sig[nc] #I Sky sigma data +real exp[nc] #I Exposure data +int bp[nc] #I Bad pixel values +int nc #I Number of columns +real cnvwt #I Sigma weight + +real a, b, s, t, rcnv, tcnv +int c, n, ncmax, nbins, bin, csky +pointer bins, skys, sigs, exps, nsky + +begin + if (line > SKB_NLSBLK(skb) * SKB_NLSPIX(skb)) + return + ncmax = min (nc, SKB_NCSBLK(skb) * SKB_NCSPIX(skb)) + + a = SKB_A(skb) + b = SKB_B(skb) + n = SKB_NCSPIX(skb) + nbins = SKB_NBINS(skb) + bins = SKB_BINS(skb) + skys = SKB_SKY(skb) + sigs = SKB_SIG(skb) + exps = SKB_EXP(skb) + nsky = SKB_NSKY(skb) + + if (IS_INDEFR(Memr[exps])) { + do c = 1, ncmax { + if (bp[c] != 0) + next + + s = sky[c] + t = sig[c] + rcnv = cnv[c] - s + tcnv = t / cnvwt + bin = a * rcnv / tcnv + b + if (bin < 1 || bin > nbins) + next + + csky = (c-1) / n + bin = bins + csky * nbins + bin - 1 + Memi[bin] = Memi[bin] + 1 + Memr[skys+csky] = Memr[skys+csky] + s + Memr[sigs+csky] = Memr[sigs+csky] + t + Memi[nsky+csky] = Memi[nsky+csky] + 1 + } + } else { + do c = 1, ncmax { + if (bp[c] != 0) + next + + s = sky[c] + t = sig[c] + rcnv = cnv[c] - s + tcnv = t / cnvwt + bin = a * rcnv / tcnv + b + if (bin < 1 || bin > nbins) + next + + csky = (c-1) / n + bin = bins + csky * nbins + bin - 1 + Memi[bin] = Memi[bin] + 1 + Memr[skys+csky] = Memr[skys+csky] + s + Memr[sigs+csky] = Memr[sigs+csky] + t + Memr[exps+csky] = Memr[exps+csky] + exp[c] + Memi[nsky+csky] = Memi[nsky+csky] + 1 + } + } + + # Evaluate histogram sky values if all lines have been accumulated. + n = mod (line, SKB_NLSPIX(skb)) + if (n == 0 || n + blkstep > SKB_NLSPIX(skb)) { + n = SKB_NCSBLK(skb) + call skb_blkeval (Memi[bins], nbins, a, b, Memr[skys], Memr[sigs], + Memr[exps], Memi[nsky], n, SKB_NSKYMIN(skb), SKB_NAV(skb), + SKB_HISTWT(skb), SKB_SIGFAC(skb)) + + # Initialize for accumulation of next line of blocks. + SKB_SKY(skb) = skys + n + SKB_SIG(skb) = sigs + n + if (!IS_INDEFR(Memr[exps])) + call aclrr (Memr[exps], n) + call aclri (Memi[nsky], n) + call aclri (Memi[bins], n*nbins) + } +end + + +# SKB_BLKEVAL -- Evaluate sky and sigma for each histogram in line of blocks. +# Set to INDEF if there are not enough pixels in the histogram. + +procedure skb_blkeval (bins, nbins, a, b, skys, sigs, exps, nsky, ncsblk, + nskymin, nav, histwt, sigfac) + +int bins[nbins,ncsblk] #I Sky subblock bins +int nbins #I Number of bins +real a, b #I Binning coefficients +real skys[ncsblk] #U Sky sum in, sky estimate out +real sigs[ncsblk] #U Sigma sum in, sigma estimate out +real exps[ncsblk] #I Exposure sum +int nsky[ncsblk] #I Number of values in bin +int ncsblk #I Number of sky pixels per subblock +int nskymin #I Minimum number of sky pixels for good sky +int nav #I Number of bins to average +int histwt #I Histogram weighting power +real sigfac #I Sigma conversion factor from mean abs dev. + +int i, j, k, l, m, n +double sky, sig, exp, x, wt, skymean, skymed, skybin, sigbin +double sum1, sum2, sum3 + +begin +# do i = 1, ncsblk { +# do j = 1, nbins { +# call printf ("%d\n") +# call pargi (bins[j,i]) +# } +# } + m = nav / 2 + do i = 1, ncsblk { + n = nsky[i] + if (n < nskymin) { + skys[i] = INDEFR + sigs[i] = INDEFR + next + } + + sky = skys[i] / n + sig = sigs[i] / n + if (!IS_INDEFR(exps[1])) { + exp = exps[i] / n + exps[i] = exp + } else + exp = 1 + + # Compute mean and median using a power weighting of the histogram. + sum1 = 0. + sum2 = 0. + sum3 = 0. + k = ncsblk + 1 + call aclri (bins[1,k], nbins) + do j = 1, nbins { + n = bins[j,i] + do l = max(1,j-m), min (nbins,j+m) + bins[l,k] = bins[l,k] + n + } + n = nsky[i] + switch (histwt) { + case 1: + do j = 1, nbins { + wt = real (bins[j,k]) / n + x = j + sum1 = sum1 + wt * x + sum2 = sum2 + wt + } + sum2 = sum2 + x = 0 + do j = 1, nbins { + wt = real (bins[j,k]) / n + sum3 = sum3 + wt + x + if (sum3 >= sum2) + break + x = wt + } + case 2: + do j = 1, nbins { + wt = real (bins[j,k]) / n + wt = wt * wt + x = j + sum1 = sum1 + wt * x + sum2 = sum2 + wt + } + sum2 = sum2 + x = 0 + do j = 1, nbins { + wt = real (bins[j,k]) / n + wt = wt * wt + sum3 = sum3 + wt + x + if (sum3 >= sum2) + break + x = wt + } + case 3: + do j = 1, nbins { + wt = real (bins[j,k]) / n + wt = wt * wt * wt + x = j + sum1 = sum1 + wt * x + sum2 = sum2 + wt + } + sum2 = sum2 + x = 0 + do j = 1, nbins { + wt = real (bins[j,k]) / n + wt = wt * wt * wt + sum3 = sum3 + wt + x + if (sum3 >= sum2) + break + x = wt + } + case 4: + do j = 1, nbins { + wt = real (bins[j,k]) / n + wt = wt * wt + wt = wt * wt + x = j + sum1 = sum1 + wt * x + sum2 = sum2 + wt + } + sum2 = sum2 + x = 0 + do j = 1, nbins { + wt = real (bins[j,k]) / n + wt = wt * wt + wt = wt * wt + sum3 = sum3 + wt + x + if (sum3 >= sum2) + break + x = wt + } + default: + do j = 1, nbins { + wt = real (bins[j,k]) / n + wt = wt ** histwt + x = j + sum1 = sum1 + wt * x + sum2 = sum2 + wt + } + sum2 = sum2 + x = 0 + do j = 1, nbins { + wt = real (bins[j,k]) / n + wt = wt ** histwt + sum3 = sum3 + wt + x + if (sum3 >= sum2) + break + x = wt + } + } + skymean = sum1 / sum2 + skymed = j - (sum3 - sum2) / (wt + x) + #skybin = skymean - max (0D0, 3 * (skymean - skymed)) + skybin = skymean - 3 * (skymean - skymed) + #skybin = skymean + skys[i] = ((skybin + 0.5 - b) / a) * sig + sky + + sum1 = 0. + sum2 = 0. + do j = 1, nbins { + wt = bins[j,k] + x = abs (j - skybin) + sum1 = sum1 + wt * x + sum2 = sum2 + wt + } + sigbin = sum1 / sum2 + sigs[i] = sigbin / a * sig * sqrt (exp) * sigfac + } +end + + +# SKB_UPDATE -- Update the sky and sigma maps using the block values. + +procedure skb_update (skb, dosky, dosig, im, skyname, signame, + skymap, sigmap, logfd) + +pointer skb #I Sky block structure +bool dosky #I Compute sky +bool dosig #I Compute sigma +pointer im #I Image pointer +char skyname[ARB] #I Output sky map name +char signame[ARB] #I Output sigma map name +pointer skymap #U Sky map pointer +pointer sigmap #U Sigma map pointer +int logfd #I Log file descriptor + +bool skydebug, sigdebug +pointer sp, fname, tmp, map_open() +errchk skb_wmap, skb_grow, skb_merge, skb_wmap, map_close, map_open + +begin + call smark (sp) + call salloc (fname, SZ_FNAME, TY_CHAR) + + + if (dosky) { + skydebug = false + if (skydebug) + call skb_wmap ("skydebug.fits", im, SKB_SKYS(skb), + SKB_NCSBLK(skb), SKB_NLSBLK(skb), SKB_NCSPIX(skb), + SKB_NLSPIX(skb), 0., NULL) + + # Grow subblocks contaminated by large objects. + call skb_grow (SKB_SKYS(skb), SKB_NCSBLK(skb), SKB_NLSBLK(skb), + SKB_GROW(skb)) + + if (skydebug) + call skb_wmap ("skydebug.fits", im, SKB_SKYS(skb), + SKB_NCSBLK(skb), SKB_NLSBLK(skb), SKB_NCSPIX(skb), + SKB_NLSPIX(skb), 0., NULL) + + # Merge sky from subblocks and interpolate missing regions. + call skb_merge (Memr[SKB_SKYS(skb)], SKB_NCSBLK(skb), + SKB_NLSBLK(skb), Memr[SKB_SKYS(skb)], SKB_NCBLK(skb), + SKB_NLBLK(skb)) + + # Write block maps and map them with the MAPIO interface. + # If no name is given then use a temporary image. + if (skyname[1] == EOS) { + call mktemp ("tmpsky", Memc[fname], SZ_FNAME) + call skb_wmap (Memc[fname], im, SKB_SKYS(skb), + SKB_NCBLK(skb), SKB_NLBLK(skb), SKB_NCPIX(skb), + SKB_NLPIX(skb), INDEFR, NULL) + } else { + call strcpy (skyname, Memc[fname], SZ_FNAME) + call skb_wmap (Memc[fname], im, SKB_SKYS(skb), + SKB_NCBLK(skb), SKB_NLBLK(skb), SKB_NCPIX(skb), + SKB_NLPIX(skb), INDEFR, logfd) + } + tmp = skymap + iferr (skymap = map_open (Memc[fname], im)) + skymap = NULL + if (skymap == NULL) { + skymap = tmp + call error (1, "Could not update sky") + } + call map_close (tmp) + if (skyname[1] == EOS) + call map_seti (skymap, "delete", YES) + } + + if (dosig) { + sigdebug = false + if (sigdebug) + call skb_wmap ("sigdebug.fits", im, SKB_SIGS(skb), + SKB_NCSBLK(skb), SKB_NLSBLK(skb), SKB_NCSPIX(skb), + SKB_NLSPIX(skb), 0., NULL) + + # Grow subblocks contaminated by large objects. + call skb_grow (SKB_SIGS(skb), SKB_NCSBLK(skb), SKB_NLSBLK(skb), + SKB_GROW(skb)) + + if (sigdebug) + call skb_wmap ("sigdebug.fits", im, SKB_SIGS(skb), + SKB_NCSBLK(skb), SKB_NLSBLK(skb), SKB_NCSPIX(skb), + SKB_NLSPIX(skb), 0., NULL) + + # Merge sky sigma from subblocks and interpolate missing regions. + call skb_merge (Memr[SKB_SIGS(skb)], SKB_NCSBLK(skb), + SKB_NLSBLK(skb), Memr[SKB_SIGS(skb)], SKB_NCBLK(skb), + SKB_NLBLK(skb)) + + # Write block maps and map them with the MAPIO interface. + # If no name is given then use a temporary image. + if (signame[1] == EOS) { + call mktemp ("tmpsig", Memc[fname], SZ_FNAME) + call skb_wmap (Memc[fname], im, SKB_SIGS(skb), + SKB_NCBLK(skb), SKB_NLBLK(skb), SKB_NCPIX(skb), + SKB_NLPIX(skb), INDEFR, NULL) + } else { + call strcpy (signame, Memc[fname], SZ_FNAME) + call skb_wmap (Memc[fname], im, SKB_SIGS(skb), + SKB_NCBLK(skb), SKB_NLBLK(skb), SKB_NCPIX(skb), + SKB_NLPIX(skb), INDEFR, logfd) + } + tmp = sigmap + iferr (sigmap = map_open (Memc[fname], im)) + sigmap = NULL + if (sigmap == NULL) { + sigmap = tmp + call error (1, "Could not update sky sigma") + } + call map_close (tmp) + if (signame[1] == EOS) + call map_seti (sigmap, "delete", YES) + } + + call sfree (sp) +end + + +# SKB_GROW -- Grow around subblocks with insufficient data. + +procedure skb_grow (sky, nc, nl, grow) + +pointer sky # Pointer to real sky array to be grown +int nc, nl # Size of sky array +real grow # Grow radius + +int i, j, k, l1, l2, ngrow, nbufs +real grow2, val1, val2, y2 +pointer buf, buf1, buf2, ptr +errchk calloc + +begin + # Initialize. + ngrow = int (grow) + grow2 = grow * grow + nbufs = min (1 + 2 * ngrow, nl) + call calloc (buf, nc*nbufs, TY_REAL) + + l1 = 1; l2 = 1 + while (l1 <= nl) { + buf1 = sky + (l1 - 1) * nc + buf2 = buf + mod (l1, nbufs) * nc + do i = 1, nc { + val1 = Memr[buf1] + val2 = Memr[buf2] + if (IS_INDEFR(val1)) { + do j = max(1,l1-ngrow), min (nl,l1+ngrow) { + ptr = buf + mod (j, nbufs) * nc - 1 + y2 = (j - l1) ** 2 + do k = max(1,i-ngrow), min (nc,i+ngrow) { + if ((k-i)**2 + y2 > grow2) + next + Memr[ptr+k] = INDEFR + } + } + } else if (!IS_INDEFR(val2)) + Memr[buf2] = val1 + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + + if (l1 > ngrow) { + while (l2 <= nl) { + buf1 = sky + (l2 - 1) * nc + buf2 = buf + mod (l2, nbufs) * nc + do i = 1, nc { + Memr[buf1] = Memr[buf2] + Memr[buf2] = 0 + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + l2 = l2 + 1 + if (l1 != nl) + break + } + } + l1 = l1 + 1 + } + + call mfree (buf, TY_REAL) +end + + +# SKB_MERGE -- Merge subblock into blocks. +# Use average of subblocks with minimum and maximum excluded. + +procedure skb_merge (in, ncin, nlin, out, ncout, nlout) + +real in[ncin,nlin] +int ncin, nlin +real out[ncout,nlout] +int ncout, nlout + +int ncs, nls +int i, i1, i2, iout, j, j1, j2, jout, n, nindef +real val, sum, minval, maxval +pointer work + +begin + # Number of input subblocks per output block. + ncs = nint (real (ncin) / ncout) + nls = nint (real (nlin) / nlout) + + nindef = 0 + j2 = 0; jout = 0 + do j1 = 1, nlin, nls { + jout = jout + 1 + j2 = min (nlin, j2 + nls) + i2 = 0; iout = 0 + do i1 = 1, ncin, ncs { + iout = iout + 1 + i2 = min (ncin, i2 + ncs) + + n = 0 + sum = 0. + minval = MAX_REAL + maxval = -MAX_REAL + do j = j1, j2 { + do i = i1, i2 { + if (IS_INDEFR(in[i,j])) + next + val = in[i,j] + sum = sum + val + minval = min (val, minval) + maxval = max (val, maxval) + n = n + 1 + } + } + if (n > 2) + out[iout,jout] = (sum - minval - maxval) / (n - 2) + else if (n >= min (ncs, nls)) + out[iout,jout] = sum / n + else { + out[iout,jout] = INDEFR + nindef = nindef + 1 + } + } + } + + # Interpolate to fill in blocks with no sky data. + if (nindef > 0) { + call malloc (work, ncout*nlout, TY_REAL) + call interp2 (out, Memr[work], ncout, nlout) + call amovr (Memr[work], out, ncout*nlout) + call mfree (work, TY_REAL) + } +end + + +## SKB_ESTIMATE -- Estimate of sky in block from subblocks. +## Use order selection. +# +#procedure skb_merge (in, ncin, nlin, out, ncout, nlout, select) +# +#real in[ncin,nlin] +#int ncin, nlin +#real out[ncout,nlout] +#int ncout, nlout +#real select # Selection fraction +# +#int ncs, nls +#int i, i1, i2, iout, j, j1, j2, jout, n, nindef, nselect +#pointer sp, work, ptr +#real asokr() +# +#begin +# # Number of input subblocks per output block. +# ncs = nint (real (ncin) / ncout) +# nls = nint (real (nlin) / nlout) +# +# call smark (sp) +# call salloc (work, ncs*nls, TY_REAL) +# +# nindef = 0 +# j2 = 0; jout = 0 +# do j1 = 1, nlin, nls { +# jout = jout + 1 +# j2 = min (nlin, j2 + nls) +# i2 = 0; iout = 0 +# do i1 = 1, ncin, ncs { +# iout = iout + 1 +# i2 = min (ncin, i2 + ncs) +# ptr = work +# do j = j1, j2 { +# do i = i1, i2 { +# if (IS_INDEFR(in[i,j])) +# next +# Memr[ptr] = in[i,j] +# ptr = ptr + 1 +# } +# } +# n = ptr - work +# if (n >= min (ncs, nls)) { +# nselect = nint (select * (n - 1)) + 1 +# out[iout,jout] = asokr (Memr[work], n, nselect) +# } else { +# out[iout,jout] = INDEFR +# nindef = nindef + 1 +# } +# } +# } +# +# # Interpolate to fill in blocks with no sky data. +# if (nindef > 0) { +# call salloc (work, ncout*nlout, TY_REAL) +# call interp2 (out, Memr[work], ncout, nlout) +# call amovr (Memr[work], out, ncout*nlout) +# } +# +# call sfree (sp) +#end + + +# SKB_WMAP -- Write map from block data. + +procedure skb_wmap (name, imref, data, ncblk, nlblk, ncpix, nlpix, blank, logfd) + +char name[ARB] #I Output name +pointer imref #I Reference image pointer +pointer data #I Block image data +int ncblk, nlblk #I Block image dimensions +int ncpix, nlpix #I Number of reference image pixels per block +real blank #I Blank value +int logfd #I Log file descriptor + +bool strne() +int i, j, imaccess(), strlen(), stridxs() +real a[2] +pointer sp, title, str +pointer im, mw, buf, immap(), impl2r(), mw_openim() +errchk immap, imrename + +begin + call smark (sp) + call salloc (title, SZ_IMTITLE, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Create title for new image or to check for updating. + call sprintf (Memc[title], SZ_IMTITLE, "Sky for ") + i = strlen (Memc[title]) + call imstats (imref, IM_IMAGENAME, Memc[title+i], SZ_IMTITLE-i) + + iferr { + im = NULL; mw = NULL + +# # Check for existing image and rename. +# if (imaccess (name, 0) == YES) { +# j = strlen (name) +# call malloc (fname, j+SZ_FNAME, TY_CHAR) +# i = strldxs (".", name) - 1 +# if (i < 0) +# i = j +# do j = 1, ARB { +# call strcpy (name, Memc[fname], i) +# call sprintf (Memc[fname+i], SZ_FNAME, "%d%s") +# call pargi (j) +# call pargstr (name[i+1]) +# if (imaccess (Memc[fname], 0) == YES) +# next +# call imrename (name, Memc[fname]) +# break +# } +# call mfree (fname, TY_CHAR) +# } + + if (imaccess (name, 0) == NO) { + if (logfd != NULL) { + call strcpy (name, Memc[str], SZ_FNAME) + i = stridxs (",", Memc[str]) + if (i > 0) { + Memc[str+i-1] = ']' + Memc[str+i] = EOS + } + call fprintf (logfd, " Write sky map: %s\n") + call pargstr (Memc[str]) + } + buf = immap (name, NEW_COPY, imref); im = buf + IM_PIXTYPE(im) = TY_REAL + IM_LEN(im,1) = ncblk + IM_LEN(im,2) = nlblk + call strcpy (Memc[title], IM_TITLE(im), SZ_IMTITLE) + iferr (call imdelf (im, "BPM")) + ; + iferr (call imdelf (im, "DATASEC")) + ; + iferr (call imdelf (im, "TRIMSEC")) + ; + + do i = 1, nlblk { + buf = impl2r(im,i) + call amovr (Memr[data+(i-1)*ncblk], Memr[buf], ncblk) + if (!IS_INDEFR(blank)) { + do j = 1, ncblk + if (IS_INDEFR(Memr[buf+j-1])) + Memr[buf+j-1] = blank + } + } + + # Update the WCS. + mw = mw_openim (imref) + a[1] = 1. / ncpix + a[2] = 1. / nlpix + call mw_scale (mw, a, 3) + a[1] = 0.5 + a[2] = 0.5 + call mw_shift (mw, a, 3) + call mw_saveim (mw, im) + } else { + if (logfd != NULL) { + call strcpy (name, Memc[str], SZ_FNAME) + i = stridxs (",", Memc[str]) + if (i > 0) { + Memc[str+i-1] = ']' + Memc[str+i] = EOS + } + call fprintf (logfd, " Update sky map: %s\n") + call pargstr (Memc[str]) + } + buf = immap (name, READ_WRITE, 0); im = buf + if (strne (IM_TITLE(im), Memc[title]) || + IM_LEN(im,1) != ncblk || IM_LEN(im,2) != nlblk) + call error (1, "Cannot update sky map") + + do i = 1, nlblk { + buf = impl2r(im,i) + call amovr (Memr[data+(i-1)*ncblk], Memr[buf], ncblk) + if (!IS_INDEFR(blank)) { + do j = 1, ncblk + if (IS_INDEFR(Memr[buf+j-1])) + Memr[buf+j-1] = blank + } + } + } + } then + call erract (EA_WARN) + + if (mw != NULL) + call mw_close (mw) + if (im != NULL) + call imunmap (im) + call sfree (sp) +end + + +# INTERP2 -- Interpolate 2D array by averaging 1D interpolations along lines +# and columns. It is an error if there is no data to interpolate. + +procedure interp2 (in, out, nc, nl) + +real in[nc,nl] # Input data +real out[nc,nl] # Output data (not the same as input) +int nc, nl # Size of data + +int i, j, k1, k2, nerr +pointer sp, flags, buf + +begin + call smark (sp) + call salloc (flags, nl, TY_INT) + call salloc (buf, nl, TY_REAL) + + call amovki (OK, Memi[flags], nl) + + # Interpolate along lines. Flag lines with no data. + nerr = 0 + do i = 1, nl + iferr (call interp1 (in[1,i], out[1,i], nc)) { + Memi[flags+i-1] = ERR + nerr = nerr + 1 + } + + if (nerr == nl) + call error (1, "No data to interpolate") + + # Interpolate along columns. Check for columns and lines with no data. + do j = 1, nc { + do i = 1, nl + Memr[buf+i-1] = in[j,i] + + ifnoerr (call interp1 (Memr[buf], Memr[buf], nl)) { + do i = 1, nl { + if (Memi[flags+i-1] == OK) + out[j,i] = (out[j,i] + Memr[buf+i-1]) / 2. + else + out[j,i] = Memr[buf+i-1] + } + } else { + do i = 1, nl { + if (Memi[flags+i-1] == ERR) { + # Find nearest line with good data. + do k1 = i-1, 1, -1 + if (Memi[flags+k1-1] == OK) + break + do k2 = i+1, nl + if (Memi[flags+k2-1] == OK) + break + if (k1 >= 1 & k2 <= nl) { + if (i - k1 < k2 - i) + out[j,i] = out[j,k1] + else + out[j,i] = out[j,k2] + } else if (k1 >= 1) + out[j,i] = out[j,k1] + else if (k2 <= nl) + out[j,i] = out[j,k2] + } + } + } + } + call sfree (sp) +end + + +# INTERP1 -- Interpolate 1D vectors. +# An error is generated if there is no data to interpolate. + +procedure interp1 (in, out, npts) + +real in[npts] # Input line +real out[npts] # Output line (may be the same as input) +int npts # Number of points in line + +int i, i1, i2, j +real v, v1, dv + +begin + i1 = 0 + i2 = 1 + do i = 1, npts { + v = in[i] + if (IS_INDEFR(v)) + next + if (i > i2) { + if (i1 > 0) { + dv = (v - v1) / (i - i1) + do j = i2, i-1 + out[j] = v + dv * (j - i) + } else { + do j = i2, i-1 + out[j] = v + } + } + out[i] = v + v1 = v + i1 = i + i2 = i1+1 + } + + if (i1 == 0) + call error (1, "No data to interpolate") + else if (i2 <= npts) { + do j = i2, npts + out[j] = v1 + } + +end |