From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001
From: Joseph Hunkeler <jhunkeler@gmail.com>
Date: Wed, 8 Jul 2015 20:46:52 -0400
Subject: Initial commit

---
 noao/imred/crutil/src/crnebula.cl | 116 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 116 insertions(+)
 create mode 100644 noao/imred/crutil/src/crnebula.cl

(limited to 'noao/imred/crutil/src/crnebula.cl')

diff --git a/noao/imred/crutil/src/crnebula.cl b/noao/imred/crutil/src/crnebula.cl
new file mode 100644
index 00000000..c5f214f9
--- /dev/null
+++ b/noao/imred/crutil/src/crnebula.cl
@@ -0,0 +1,116 @@
+# CRNEBULA -- Cosmic ray cleaning for images with fine nebular structure.
+
+procedure crnebula (input, output)
+
+file    input           {prompt="Input image"}
+file	output		{prompt="Output image"}
+file    crmask          {prompt="Cosmic ray mask"}
+file	residual	{prompt="Residual median image"}
+file	rmedresid	{prompt="Residual ring median image"}
+real    var0 = 1.       {prompt="Variance coefficient for DN^0 term", min=0.}
+real	var1 = 0.	{prompt="Variance coefficient for DN^1 term", min=0.}
+real	var2 = 0.	{prompt="Variance coefficient for DN^2 term", min=0.}
+real    sigmed = 3      {prompt="Sigma clip factor for residual median"}
+real    sigrmed = 3     {prompt="Sigma clip factor for residual ring median"}
+
+int     mbox = 5        {prompt="Median box size"}
+real    rin = 1.5       {prompt="Inner radius for ring median"}
+real    rout = 6        {prompt="Outer radius for ring median"}
+bool    verbose = no    {prompt="Verbose"}
+
+begin
+        file    in, out, med, sig, resid, rmed
+        struct  expr
+
+        # Query once for query parameters.
+        in = input
+        out = output
+
+        # Check on output images.
+	if (out != "" && imaccess (out))
+            error (1, "Output image already exists ("//out//")")
+        if (crmask != "" && imaccess (crmask))
+            error (1, "Output mask already exists ("//crmask//")")
+	if (residual != "" && imaccess (residual))
+            error (1, "Output residual image already exists ("//residual//")")
+	if (rmedresid != "" && imaccess (rmedresid))
+            error (1,
+	       "Output ring median difference already exists ("//rmedresid//")")
+
+        # Create median results.
+	med = mktemp ("cr")
+	sig = mktemp ("cr")
+	resid = residual
+	if (resid == "")
+	    resid = mktemp ("cr")
+	if (verbose)
+	    printf ("Creating CRMEDIAN results\n")
+	crmedian (in, "", crmask="", median=med, sigma=sig, residual=resid,
+	    var0=var0, var1=var1, var2=var2, lsigma=100., hsigma=sigmed,
+	    ncmed=mbox, nlmed=mbox, ncsig=25, nlsig=25)
+
+        # Create ring median filtered image.
+        rmed = mktemp ("cr")
+	rmedian (in, rmed, rin, rout, ratio=1., theta=0., zloreject=INDEF,
+	    zhireject=INDEF, boundary="wrap", constant=0., verbose=verbose)
+
+	# Create output images.
+	if (rmedresid != "") {
+	    printf ("(a-b)/c\n") | scan (expr)
+	    imexpr (expr, rmedresid, med, rmed, sig, dims="auto",
+		intype="auto", outtype="real", refim="auto", bwidth=0,
+		btype="nearest", bpixval=0., rangecheck=yes, verbose=no,
+		exprdb="none")
+	    imdelete (rmed, verify-)
+	    imdelete (sig, verify-)
+	    if (out != "") {
+		if (verbose)
+		    printf ("Create output image %s\n", out)
+		printf ("((a<%.3g)||(abs(b)>%.3g)) ? c : d\n",
+		    sigmed, sigrmed) | scan (expr)
+		imexpr (expr, out, resid, rmedresid, in, med, dims="auto",
+		    intype="auto", outtype="auto", refim="auto", bwidth=0,
+		    btype="nearest", bpixval=0., rangecheck=yes, verbose=no,
+		    exprdb="none")
+	    }
+	    if (crmask != "") {
+		if (verbose)
+		    printf ("Create cosmic ray mask %s\n", crmask)
+		printf ("((a<%.3g)||(abs(b)>%.3g)) ? 0 : 1\n",
+		    sigmed, sigrmed) | scan (expr)
+		set imtype = "pl"
+		imexpr (expr, crmask, resid, rmedresid, dims="auto",
+		    intype="auto", outtype="short", refim="auto", bwidth=0,
+		    btype="nearest", bpixval=0., rangecheck=yes, verbose=no,
+		    exprdb="none")
+	    }
+	    imdelete (med, verify-)
+	} else {
+	    if (out != "") {
+		if (verbose)
+		    printf ("Create output image %s\n", out)
+		printf ("((a<%.3g)||(abs((b-c)/d)>%.3g)) ? e : b\n",
+		    sigmed, sigrmed) | scan (expr)
+		imexpr (expr, out, resid, med, rmed, sig, in, dims="auto",
+		    intype="auto", outtype="auto", refim="auto", bwidth=0,
+		    btype="nearest", bpixval=0., rangecheck=yes, verbose=no,
+		    exprdb="none")
+	    }
+	    if (crmask != "") {
+		if (verbose)
+		    printf ("Create cosmic ray mask %s\n", crmask)
+		printf ("((a<%.3g)||(abs((b-c)/d)>%.3g)) ? 0 : 1\n",
+		    sigmed, sigrmed) | scan (expr)
+		set imtype = "pl"
+		imexpr (expr, crmask, resid, med, rmed, sig, dims="auto",
+		    intype="auto", outtype="short", refim="auto", bwidth=0,
+		    btype="nearest", bpixval=0., rangecheck=yes, verbose=no,
+		    exprdb="none")
+	    }
+	    imdelete (med, verify-)
+	    imdelete (sig, verify-)
+	    imdelete (rmed, verify-)
+	}
+	if (residual == "")
+	    imdelete (resid, verify-)
+end
-- 
cgit