Repatch (from linux) of OSX IRAF

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