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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
|
include <imhdr.h>
include "rskysub.h"
# T_RSKYSUB -- Sky subtract a set of input images using image scaling and
# a running statistics compution
procedure t_rskysub()
pointer sp, imasks, str
pointer rs
int inlist, imsklist, outlist, omsklist, hmsklist, sclist, tmplist
bool msk_invert, useimasks, cache, verbose
real clgetr()
int imtopenp(), imtopen(), imtlen(), fntopnb(), fntlenb()
int clgeti(), btoi(), clgwrd(), rs_imlist(), rs_olist(), rs_omlist()
bool clgetb(), strne()
begin
# Allocate working space.
call smark (sp)
call salloc (imasks, SZ_FNAME, TY_CHAR)
call salloc (str, SZ_FNAME, TY_CHAR)
# Open the input image list. Make this a test versus nmin ?
inlist = imtopenp ("input")
if (imtlen (inlist) <= 0) {
call eprintf ("The input image list is empty\n")
call imtclose (inlist)
call sfree (sp)
return
}
# Open the output image list. The number of output images must be
# zero equal to the number of input images.
call clgstr ("output", Memc[str], SZ_FNAME)
outlist = rs_olist (inlist, Memc[str], "default", "sub")
if (imtlen (outlist) > 0 && imtlen(outlist) != imtlen(inlist)) {
call eprintf ("The output mask and image lists don't match\n")
call imtclose (outlist)
call imtclose (inlist)
call sfree (sp)
return
}
# Open the input mask list.
call clgstr ("imasks", Memc[imasks], SZ_FNAME)
if (Memc[imasks] == '^') {
#imsklist = imtopen (Memc[imasks+1])
imsklist = rs_imlist (inlist, Memc[imasks+1], "default", "obm")
msk_invert = true
} else {
#imsklist = imtopen (Memc[imasks])
imsklist = rs_imlist (inlist, Memc[imasks], "default", "obm")
msk_invert = false
}
if (imtlen (imsklist) > 1 && imtlen (imsklist) != imtlen (inlist)) {
call eprintf ("The input mask and image lists don't match\n")
call imtclose (imsklist)
call imtclose (outlist)
call imtclose (inlist)
call sfree (sp)
return
}
# Open the output mask list. The number of output masks must be
# zero equal to the number of input images.
call clgstr ("omasks", Memc[str], SZ_FNAME)
omsklist = rs_omlist (inlist, Memc[str], "default", "skm")
if (imtlen (omsklist) > 0 && imtlen(omsklist) != imtlen(inlist)) {
call eprintf ("The output mask and image lists don't match\n")
call imtclose (omsklist)
call imtclose (imsklist)
call imtclose (outlist)
call imtclose (inlist)
call sfree (sp)
return
}
# Open the output holes mask list. The number of output holes masks
# must be zero equal to the number of input images.
call clgstr ("hmasks", Memc[str], SZ_FNAME)
hmsklist = rs_omlist (inlist, Memc[str], "default", "hom")
if (imtlen (hmsklist) > 0 && imtlen(hmsklist) != imtlen(inlist)) {
call eprintf ("The holes mask and image lists don't match\n")
call imtclose (hmsklist)
call imtclose (omsklist)
call imtclose (imsklist)
call imtclose (outlist)
call imtclose (inlist)
call sfree (sp)
return
}
# Allocate the sky subtraction structure
call malloc (rs, LEN_RSKYSUB, TY_STRUCT)
# Get the scaling factor computation method.
RS_RESCALE(rs) = btoi(clgetb ("rescale"))
call clgstr ("scale", RS_ISCALES(rs), SZ_FNAME)
sclist = fntopnb (RS_ISCALES(rs), NO)
if (fntlenb (sclist) > 1 && fntlenb (sclist) != imtlen (inlist)) {
call eprintf ("The scaling factor and image lists don't match\n")
call fntclsb (sclist)
call imtclose (hmsklist)
call imtclose (omsklist)
call imtclose (imsklist)
call imtclose (outlist)
call imtclose (inlist)
call sfree (sp)
return
}
# If the scaling algorith is not "median" then new output masks
# cannot be created.
if (strne (RS_ISCALES(rs), "median")) {
call imtclose (omsklist)
omsklist = imtopen ("")
}
call clgstr ("skyscale", RS_KYFSCALE(rs), SZ_FNAME)
# Get statisitics computation parameters.
useimasks = clgetb ("useimasks")
call clgstr ("statsec", RS_STATSEC(rs), SZ_FNAME)
RS_LOWER(rs) = clgetr ("lower")
RS_UPPER(rs) = clgetr ("upper")
RS_MAXITER(rs) = clgeti ("maxiter")
RS_LNSIGREJ(rs) = clgetr ("lnsigrej")
RS_UNSIGREJ(rs) = clgetr ("unsigrej")
RS_BINWIDTH(rs) = clgetr ("binwidth")
if (RS_MAXITER(rs) > 0 && IS_INDEFR(RS_LNSIGREJ(rs)) &&
IS_INDEFR(RS_UNSIGREJ(rs)))
RS_MAXITER(rs) = 0
# Get the sky subtraction parameters
RS_RESUBTRACT(rs) = btoi(clgetb ("resubtract"))
RS_COMBINE(rs) = clgwrd ("combine", Memc[str], SZ_FNAME, RS_COMBINESTR)
RS_NCOMBINE(rs) = clgeti ("ncombine")
RS_NMIN(rs) = clgeti ("nmin")
if (RS_NMIN(rs) <= 0 || RS_NMIN(rs) > RS_NCOMBINE(rs)) {
RS_NMIN(rs) = RS_NCOMBINE(rs)
call eprintf ("Warning: resetting nmin to %d\n")
call pargi (RS_NMIN(rs))
}
# Get starting values for the rejection parameters. These may have
# to be adjusted if image masking is enabled and for cases where
# the number of combined images is greater then equal to nmin but
# less than ncombine.
RS_NLOREJ(rs) = clgeti ("nlorej")
RS_NHIREJ(rs) = clgeti ("nhirej")
switch (RS_COMBINE(rs)) {
case RS_MEAN:
if ((RS_NMIN(rs) - RS_NLOREJ(rs) - RS_NHIREJ(rs)) < 1) {
call eprintf ("Too many rejected pixels\n")
call fntclsb (sclist)
call imtclose (hmsklist)
call imtclose (omsklist)
call imtclose (imsklist)
call imtclose (outlist)
call imtclose (inlist)
call sfree (sp)
return
}
case RS_MEDIAN:
if (mod (RS_NCOMBINE(rs), 2) == 0) {
RS_NLOREJ(rs) = RS_NCOMBINE(rs) / 2 - 1
RS_NHIREJ(rs) = RS_NCOMBINE(rs) / 2 - 1
} else {
RS_NLOREJ(rs) = RS_NCOMBINE(rs) / 2
RS_NHIREJ(rs) = RS_NCOMBINE(rs) / 2
}
default:
}
RS_BLANK(rs) = clgetr ("blank")
call clgstr ("skysub", RS_KYSKYSUB(rs), SZ_FNAME)
call clgstr ("holes", RS_KYHMASK(rs), SZ_FNAME)
cache = clgetb ("cache")
verbose = clgetb ("verbose")
# Compute the sky statistics and optionally the output sky masks.
if (useimasks) {
call rs_stats (inlist, imsklist, omsklist, sclist, rs, msk_invert,
cache, verbose)
} else {
tmplist = imtopen ("")
call rs_stats (inlist, tmplist, omsklist, sclist, rs, msk_invert,
cache, verbose)
call imtclose (tmplist)
}
# Do the sky subtraction with or without image masking and with or
# without bad pixel rejection. Unmasked image medians can be handled
# by setting the high and low pixel rejection parameters appropriately.
# Masked image means and medians may require dynaimc altering of the
# high and low rejection parameters.
switch (RS_COMBINE(rs)) {
case RS_MEAN, RS_MEDIAN:
if (imtlen (omsklist) > 0) {
if (RS_NLOREJ(rs) > 0 || RS_NHIREJ(rs) > 0)
# Choose which of the two routines to use later based
# on timing tests.
#call rs_prmsub (inlist, omsklist, outlist, hmsklist, rs,
#msk_invert, cache, verbose)
call rs_prrmsub (inlist, omsklist, outlist, hmsklist, rs,
msk_invert, cache, verbose)
else
call rs_pmsub (inlist, omsklist, outlist, hmsklist, rs,
msk_invert, cache, verbose)
} else if (imtlen (imsklist) > 0) {
if (RS_NLOREJ(rs) > 0 || RS_NHIREJ(rs) > 0)
# Choose which of the two routines to use later based on
# timing tests.
#call rs_prmsub (inlist, imsklist, outlist, hmsklist, rs,
#msk_invert, cache, verbose)
call rs_prrmsub (inlist, imsklist, outlist, hmsklist, rs,
msk_invert, cache, verbose)
else
call rs_pmsub (inlist, imsklist, outlist, hmsklist, rs,
msk_invert, cache, verbose)
} else {
if (RS_NLOREJ(rs) > 0 || RS_NHIREJ(rs) > 0)
# Choose which of the two routines to use later based on
# timing tests.
#call rs_rmsub (inlist, outlist, rs, cache, verbose)
call rs_rrmsub (inlist, outlist, rs, cache, verbose)
else
call rs_msub (inlist, outlist, rs, cache, verbose)
}
default:
;
}
# Close image and file lists.
call fntclsb (sclist)
call imtclose (hmsklist)
call imtclose (imsklist)
call imtclose (omsklist)
call imtclose (outlist)
call imtclose (inlist)
call mfree (rs, TY_STRUCT)
call sfree (sp)
end
|
