1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
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
|
