Initial commit

1 files changed, 329 insertions, 0 deletions

diff --git a/noao/imred/crutil/src/t_cosmicrays.x b/noao/imred/crutil/src/t_cosmicrays.x
new file mode 100644
index 00000000..fa68a39d
--- /dev/null
+++ b/noao/imred/crutil/src/t_cosmicrays.x
@@ -0,0 +1,329 @@
+include	<error.h>
+include <imhdr.h>
+include <imset.h>
+include	<math/gsurfit.h>
+include	<gset.h>
+include	<pkg/gtools.h>
+include	"crlist.h"
+
+# T_COSMICRAYS -- Detect and remove cosmic rays in images.
+# A list of images is examined for cosmic rays which are then replaced
+# by values from neighboring pixels.  The output image may be the same
+# as the input image.  This is the top level procedure which manages
+# the input and output image data.  The actual algorithm for detecting
+# cosmic rays is in CR_FIND.
+
+procedure t_cosmicrays ()
+
+int	list1			# List of input images to be cleaned
+int	list2			# List of output images
+int	list3			# List of output bad pixel files
+real	threshold		# Detection threshold
+real	fluxratio		# Luminosity boundary for stars
+int	npasses			# Number of cleaning passes
+int	szwin			# Size of detection window
+bool	train			# Use training objects?
+pointer	savefile		# Save file for training objects
+bool	interactive		# Examine cosmic ray parameters?
+char	ans			# Answer to interactive query
+
+int	nc, nl, c, c1, c2, l, l1, l2, szhwin, szwin2
+int	i, j, k, m, ncr, ncrlast, nreplaced, flag
+pointer	sp, input, output, badpix, str, gp, gt, im, in, out
+pointer	x, y, z, w, sf1, sf2, cr, data, ptr
+
+bool	clgetb(), streq(), strne()
+char	clgetc()
+int	imtopenp(), imtlen(), imtgetim(), clpopnu(), clgfil(), clgeti()
+real	clgetr()
+pointer	immap(), impl2r(), imgs2r(), gopen(), gt_init()
+errchk	immap, impl2r, imgs2r
+errchk	cr_find, cr_examine, cr_replace, cr_plot, cr_crmask
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (badpix, SZ_FNAME, TY_CHAR)
+	call salloc (savefile, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Get the task parameters.  Check that the number of output images
+	# is either zero, in which case the cosmic rays will be removed
+	# in place, or equal to the number of input images.
+
+	list1 = imtopenp ("input")
+	list2 = imtopenp ("output")
+	i = imtlen (list1)
+	j = imtlen (list2)
+	if (j > 0 && j != i)
+	    call error (0, "Input and output image lists do not match")
+
+	list3 = clpopnu ("crmasks")
+	threshold = clgetr ("threshold")
+	fluxratio = clgetr ("fluxratio")
+	npasses = clgeti ("npasses")
+	szwin = clgeti ("window")
+	train = clgetb ("train")
+	call clgstr ("savefile", Memc[savefile], SZ_FNAME)
+	interactive = clgetb ("interactive")
+	call clpstr ("answer", "yes")
+	ans = 'y'
+
+	# Set up the graphics.
+	call clgstr ("graphics", Memc[str], SZ_LINE)
+	if (interactive) {
+	    gp = gopen (Memc[str], NEW_FILE+AW_DEFER, STDGRAPH)
+	    gt = gt_init()
+	    call gt_sets (gt, GTTYPE, "mark")
+	    call gt_sets (gt, GTXTRAN, "log")
+	    call gt_setr (gt, GTXMIN, 10.)
+	    call gt_setr (gt, GTYMIN, 0.)
+	    call gt_sets (gt, GTTITLE, "Parameters of cosmic rays candidates")
+	    call gt_sets (gt, GTXLABEL, "Flux")
+	    call gt_sets (gt, GTYLABEL, "Flux Ratio")
+	}
+
+	# Set up surface fitting.  The background points are placed together
+	# at the beginning of the arrays.  There are two surface pointers,
+	# one for using the fast refit if there are no points excluded and
+	# one for doing a full fit with points excluded.
+
+	szhwin = szwin / 2
+	szwin2 = szwin * szwin
+	call salloc (data, szwin, TY_INT)
+	call salloc (x, szwin2, TY_REAL)
+	call salloc (y, szwin2, TY_REAL)
+	call salloc (z, szwin2, TY_REAL)
+	call salloc (w, szwin2, TY_REAL)
+
+	k = 0
+	do i = 1, szwin {
+	    Memr[x+k] = i
+	    Memr[y+k] = 1
+	    k = k + 1
+	}
+	do i = 2, szwin {
+	    Memr[x+k] = szwin
+	    Memr[y+k] = i
+	    k = k + 1
+	}
+	do i = szwin-1, 1, -1 {
+	    Memr[x+k] = i
+	    Memr[y+k] = szwin
+	    k = k + 1
+	}
+	do i = szwin-1, 2, -1 {
+	    Memr[x+k] = 1
+	    Memr[y+k] = i
+	    k = k + 1
+	}
+	do i = 2, szwin-1 {
+	    do j = 2, szwin-1 {
+	        Memr[x+k] = j
+	        Memr[y+k] = i
+	        k = k + 1
+	    }
+	}
+	call aclrr (Memr[z], szwin2)
+	call amovkr (1., Memr[w], 4*szwin-4)
+	call gsinit (sf1, GS_POLYNOMIAL, 2, 2, NO, 1., real(szwin),
+	    1., real(szwin))
+	call gsinit (sf2, GS_POLYNOMIAL, 2, 2, NO, 1., real(szwin),
+	    1., real(szwin))
+	call gsfit (sf1, Memr[x], Memr[y], Memr[z], Memr[w], 4*szwin-4,
+	    WTS_USER, j)
+
+	# Process each input image.  Either work in place or create a
+	# new output image.  If an error mapping the images occurs
+	# issue a warning and go on to the next input image.
+
+	while (imtgetim (list1, Memc[input], SZ_FNAME) != EOF) {
+	    if (imtgetim (list2, Memc[output], SZ_FNAME) == EOF)
+		call strcpy (Memc[input], Memc[output], SZ_FNAME)
+	    if (clgfil (list3, Memc[badpix], SZ_FNAME) == EOF)
+		Memc[badpix] = EOS
+
+	    iferr {
+		in = NULL
+		out = NULL
+		cr = NULL
+
+		# Map the input image.  If the output image is
+		# the same as the input image work in place.
+		# Initialize IMIO to use a scrolling buffer of lines.
+
+	        if (streq (Memc[input], Memc[output])) {
+	            im = immap (Memc[input], READ_WRITE, 0)
+		} else
+		    im = immap (Memc[input], READ_ONLY, 0)
+		in = im
+
+	        nc = IM_LEN(in,1)
+	        nl = IM_LEN(in,2)
+	        if ((nl < szwin) || (nc < szwin))
+	            call error (0, "Image size is too small")
+	        call imseti (in, IM_NBUFS, szwin)
+	        call imseti (in, IM_TYBNDRY, BT_NEAREST)
+	        call imseti (in, IM_NBNDRYPIX, szhwin)
+
+		# Open the output image if needed.
+	        if (strne (Memc[input], Memc[output]))
+		    im = immap (Memc[output], NEW_COPY, in)
+		out = im
+
+		# Open a cosmic ray list structure.
+	        call cr_open (cr)
+		ncrlast = 0
+		nreplaced = 0
+
+		# Now proceed through the image line by line, scrolling
+		# the line buffers at each step.  If creating a new image
+		# also write out each line as it is read.  A procedure is
+		# called to find the cosmic ray candidates in the line
+		# and add them to the list maintained by CRLIST.
+		# Note that cosmic rays are not replaced at this point
+		# in order to allow the user to modify the criteria for
+		# a cosmic ray and review the results.
+
+	        c1 = 1-szhwin
+	        c2 = nc+szhwin
+		do i = 1, szwin-1 
+		    Memi[data+i] =
+			imgs2r (in, c1, c2, i-szhwin, i-szhwin)
+
+	        do l = 1, nl {
+		    do i = 1, szwin-1
+			Memi[data+i-1] = Memi[data+i]
+		    Memi[data+szwin-1] =
+			imgs2r (in, c1, c2, l+szhwin, l+szhwin)
+		    if (out != in)
+		        call amovr (Memr[Memi[data+szhwin]+szhwin],
+			    Memr[impl2r(out,l)], nc)
+
+	            call cr_find (cr, threshold, Memi[data],
+		        c2-c1+1, szwin, c1, l,
+		        sf1, sf2, Memr[x], Memr[y], Memr[z], Memr[w])
+	        }
+		if (interactive && train) {
+		    call cr_train (cr, gp, gt, in, fluxratio, Memc[savefile])
+		    train = false
+		}
+	        call cr_flags (cr, fluxratio)
+
+		# If desired examine the cosmic ray list interactively.
+	        if (interactive && ans != 'N') {
+		    if (ans != 'Y') {
+		        call eprintf ("%s - ")
+		            call pargstr (Memc[input])
+		        call flush (STDERR)
+		        ans = clgetc ("answer")
+		    }
+		    if ((ans == 'Y') || (ans == 'y'))
+	                call cr_examine (cr, gp, gt, in, fluxratio, 'r')
+	        }
+
+		# Now replace the selected cosmic rays in the output image.
+
+		call imflush (out)
+		call imseti (out, IM_ADVICE, RANDOM)
+	        call cr_replace (cr, ncrlast, out, nreplaced)
+
+		# Do additional passes through the data.  We work in place
+		# in the output image.  Note that we only have to look in
+		# the vicinity of replaced cosmic rays for secondary
+		# events since we've already looked at every pixel once.
+		# Instead of scrolling through the image we will extract
+		# subrasters around each replaced cosmic ray.  However,
+		# we use pointers into the subraster to maintain the same
+		# format expected by CR_FIND.
+
+	        if (npasses > 1) {
+	            if (out != in)
+	                call imunmap (out)
+	            call imunmap (in)
+		    im = immap (Memc[output], READ_WRITE, 0)
+		    in = im
+		    out = im
+	            call imseti (in, IM_TYBNDRY, BT_NEAREST)
+	            call imseti (in, IM_NBNDRYPIX, szhwin)
+
+	            for (i=2; i<=npasses; i=i+1) {
+			# Loop through each cosmic ray in the previous pass.
+		        ncr = CR_NCR(cr)
+		        do j = ncrlast+1, ncr {
+			    flag = Memi[CR_FLAG(cr)+j-1]
+			    if (flag==NO || flag==ALWAYSNO)
+			        next
+			    c = Memr[CR_COL(cr)+j-1]
+			    l = Memr[CR_LINE(cr)+j-1]
+			    c1 = max (1-szhwin, c - (szwin-1))
+			    c2 = min (nc+szhwin, c + (szwin-1))
+			    k = c2 - c1 + 1
+			    l1 = max (1-szhwin, l - (szwin-1))
+			    l2 = min (nl+szhwin, l + (szwin-1))
+
+			    # Set the line pointers off an image section
+			    # centered on a previously replaced cosmic ray.
+
+			    ptr = imgs2r (in, c1, c2, l1, l2) - k
+
+			    l1 = max (1, l - szhwin)
+			    l2 = min (nl, l + szhwin)
+			    do l = l1, l2 {
+				do m = 1, szwin
+				    Memi[data+m-1] = ptr + m * k
+				ptr = ptr + k
+
+	                        call cr_find ( cr, threshold, Memi[data],
+				    k, szwin, c1, l, sf1, sf2,
+				    Memr[x], Memr[y], Memr[z], Memr[w])
+			    }
+	                }
+		        call cr_flags (cr, fluxratio)
+
+			# Replace any new cosmic rays found.
+	                call cr_replace (cr, ncr, in, nreplaced)
+			ncrlast = ncr
+		    }
+	        }
+
+		# Output header log, log, plot, and bad pixels.
+		call sprintf (Memc[str], SZ_LINE,
+		    "Threshold=%5.1f, fluxratio=%6.2f, removed=%d")
+		    call pargr (threshold)
+		    call pargr (fluxratio)
+		    call pargi (nreplaced)
+		call imastr (out, "crcor", Memc[str])
+
+		call cr_plot (cr, in, fluxratio)
+		call cr_crmask (cr, Memc[badpix], in)
+
+	        call cr_close (cr)
+	        if (out != in)
+	            call imunmap (out)
+	        call imunmap (in)
+	    } then {
+		# In case of error clean up and go on to the next image.
+		if (in != NULL) {
+		    if (out != NULL && out != in)
+			call imunmap (out)
+		    call imunmap (in)
+		}
+		if (cr != NULL)
+		    call cr_close (cr)
+		call erract (EA_WARN)
+	    }
+	}
+
+	if (interactive) {
+	    call gt_free (gt)
+	    call gclose (gp)
+	}
+	call imtclose (list1)
+	call imtclose (list2)
+	call clpcls (list3)
+	call gsfree (sf1)
+	call gsfree (sf2)
+	call sfree (sp)
+end