diff options
Diffstat (limited to 'noao/imred/crutil/src/t_crgrow.x')
| -rw-r--r-- | noao/imred/crutil/src/t_crgrow.x | 182 |
1 files changed, 182 insertions, 0 deletions
diff --git a/noao/imred/crutil/src/t_crgrow.x b/noao/imred/crutil/src/t_crgrow.x new file mode 100644 index 00000000..7316cc76 --- /dev/null +++ b/noao/imred/crutil/src/t_crgrow.x @@ -0,0 +1,182 @@ +include <error.h> +include <imhdr.h> + +# T_CRGROW -- Grow cosmic ray mask identifications. + +procedure t_crgrow () + +int input # Input masks +int output # Output masks +real radius # Radius +int inval # Input mask value to grow +int outval # Output grown mask value + +pointer sp, inmask, outmask, temp1, temp2 +pointer im, ptr + +int imtopenp(), imtlen(), imtgetim() +bool strne() +int clgeti() +real clgetr() +pointer immap() +errchk immap, crgrow + +begin + call smark (sp) + call salloc (inmask, SZ_FNAME, TY_CHAR) + call salloc (outmask, SZ_FNAME, TY_CHAR) + call salloc (temp1, SZ_FNAME, TY_CHAR) + call salloc (temp2, SZ_FNAME, TY_CHAR) + + # Task parameters. + input = imtopenp ("input") + output = imtopenp ("output") + radius = max (0., clgetr ("radius")) + inval = clgeti ("inval") + outval = clgeti ("outval") + + if (imtlen (output) != imtlen (input)) + call error (1, "Input and output lists do not match") + + # Grow the cosmic ray masks. + while (imtgetim (input, Memc[inmask], SZ_FNAME) != EOF) { + call strcpy (Memc[inmask], Memc[outmask], SZ_FNAME) + if (imtgetim (output, Memc[outmask], SZ_FNAME) == EOF) + call error (1, "Output list ended prematurely") + if (strne (Memc[inmask], Memc[outmask])) { + call imgcluster (Memc[inmask], Memc[temp1], SZ_FNAME) + call imgcluster (Memc[outmask], Memc[temp2], SZ_FNAME) + iferr (call imcopy (Memc[temp1], Memc[temp2])) { + call erract (EA_WARN) + next + } + im = immap (Memc[inmask], READ_ONLY, TY_CHAR) + iferr (call imgstr (im, "extname", Memc[temp1], SZ_FNAME)) + call strcpy ("pl", Memc[temp1], SZ_FNAME) + call imunmap (im) + } + call xt_maskname (Memc[outmask], Memc[temp1], 0, Memc[outmask], + SZ_FNAME) + + if (radius < 1.) + next + + iferr { + im = NULL + ptr = immap (Memc[outmask], READ_WRITE, 0); im = ptr + call crgrow (im, radius, inval, outval) + } then { + call erract (EA_WARN) + if (strne (Memc[inmask], Memc[outmask])) { + if (im != NULL) { + call imunmap (im) + iferr (call imdelete (Memc[outmask])) + call erract (EA_WARN) + } + } + } + + if (im != NULL) + call imunmap (im) + } + + call imtclose (output) + call imtclose (input) + call sfree (sp) +end + + +# CRGROW -- Grow cosmic rays. + +procedure crgrow (im, grow, inval, outval) + +pointer im # Mask pointer (Read/Write) +real grow # Radius (pixels) +int inval # Input mask value for pixels to grow +int outval # Output mask value for grown pixels + +int i, j, k, l, nc, nl, ngrow, nbufs, val1, val2 +long v1[2], v2[2] +real grow2, y2 +pointer , buf, buf1, buf2, ptr + +int imgnli(), impnli() +errchk calloc, imgnli, impnli + +begin + if (grow < 1. || inval == 0) + return + + grow2 = grow * grow + ngrow = int (grow) + buf = NULL + + iferr { + if (IM_NDIM(im) > 2) + call error (1, + "Only one or two dimensional masks are allowed") + + nc = IM_LEN(im, 1) + if (IM_NDIM(im) > 1) + nl = IM_LEN(im,2) + else + nl = 1 + + # Initialize buffering. + nbufs = min (1 + 2 * ngrow, nl) + call calloc (buf, nc*nbufs, TY_INT) + + call amovkl (long(1), v1, IM_NDIM(im)) + call amovkl (long(1), v2, IM_NDIM(im)) + while (imgnli (im, buf1, v1) != EOF) { + j = v1[2] - 1 + buf2 = buf + mod (j, nbufs) * nc + do i = 1, nc { + val1 = Memi[buf1] + val2 = Memi[buf2] + if ((IS_INDEFI(inval) && val1 != 0) || val1 == inval) { + do k = max(1,j-ngrow), min (nl,j+ngrow) { + ptr = buf + mod (k, nbufs) * nc - 1 + y2 = (k - j) ** 2 + do l = max(1,i-ngrow), min (nc,i+ngrow) { + if ((l-i)**2 + y2 > grow2) + next + Memi[ptr+l] = -val1 + } + } + } else { + if (val2 >= 0) + Memi[buf2] = val1 + } + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + + if (j > ngrow) { + while (impnli (im, buf2, v2) != EOF) { + k = v2[2] - 1 + buf1 = buf + mod (k, nbufs) * nc + do i = 1, nc { + val1 = Memi[buf1] + if (val1 < 0) { + if (IS_INDEFI(outval)) + Memi[buf2] = -val1 + else + Memi[buf2] = outval + } else + Memi[buf2] = val1 + Memi[buf1] = 0. + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + if (j != nl) + break + } + } + } + } then + call erract (EA_ERROR) + + if (buf != NULL) + call mfree (buf, TY_INT) +end |