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
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <imhdr.h>
# IMSUM -- Sum or average images with optional high and low pixel rejection.
procedure t_imsum ()
int list # Input image list
pointer image # Output image
pointer hparams # Header parameter list
pointer option # Output option
int pixtype # Output pixel datatype
int calctype # Internal calculation type
real low_reject # Number or frac of low pix to reject
real high_reject # Number or frac of high pix to reject
int i, nimages, nlow, nhigh
pointer sp, str, im_in, im_out
bool clgetb(), streq()
real clgetr()
int imtopenp(), imtlen(), imtgetim(), clgwrd()
pointer immap()
errchk imsum_set, immap, imunmap
begin
# Get the input image list. Check that there is at least 1 image.
list = imtopenp ("input")
nimages = imtlen (list)
if (nimages < 1) {
call imtclose (list)
call error (0, "No input images in list")
}
# Allocate strings and get the parameters.
call smark (sp)
call salloc (image, SZ_FNAME, TY_CHAR)
call salloc (hparams, SZ_LINE, TY_CHAR)
call salloc (option, SZ_LINE, TY_CHAR)
i = clgwrd ("option", Memc[option], SZ_LINE, "|sum|average|median|")
if (streq (Memc[option], "median")) {
nlow = nimages / 2
nhigh = nimages - nlow - 1
} else {
# If the rejection value is less than 1 then it is a fraction of the
# input images otherwise it is the number of pixels to be rejected.
low_reject = clgetr ("low_reject")
high_reject = clgetr ("high_reject")
if (low_reject < 1.)
nlow = low_reject * nimages
else
nlow = low_reject
if (high_reject < 1.)
nhigh = high_reject * nimages
else
nhigh = high_reject
if (nlow + nhigh >= nimages) {
call sfree (sp)
call imtclose (list)
call error (0, "Number of pixels rejected >= number of images")
}
}
call clgstr ("hparams", Memc[hparams], SZ_LINE)
# Map the output image and set the title and pixel type.
# Check all images have the same number and length of dimensions.
call imsum_set (list, pixtype, calctype)
i = imtgetim (list, Memc[image], SZ_FNAME)
im_in = immap (Memc[image], READ_ONLY, 0)
call clgstr ("output", Memc[image], SZ_FNAME)
im_out = immap (Memc[image], NEW_COPY, im_in)
call new_title ("title", im_out)
IM_PIXTYPE (im_out) = pixtype
call imtrew (list)
# Print verbose info.
if (clgetb ("verbose")) {
call salloc (str, SZ_LINE, TY_CHAR)
call printf ("IMSUM:\n")
call printf (" Input images:\n")
while (imtgetim (list, Memc[str], SZ_LINE) != EOF) {
call printf (" %s\n")
call pargstr (Memc[str])
}
call imtrew (list)
call printf (" Output image: %s\n")
call pargstr (Memc[image])
call printf (" Header parameters: %s\n")
call pargstr (Memc[hparams])
call printf (" Output pixel datatype: %s\n")
call dtstring (pixtype, Memc[str], SZ_FNAME)
call pargstr (Memc[str])
call printf (" Calculation type: %s\n")
call dtstring (calctype, Memc[str], SZ_FNAME)
call pargstr (Memc[str])
call printf (" Option: %s\n")
call pargstr (Memc[option])
call printf (" Low rejection: %d\n High rejection: %d\n")
call pargi (nlow)
call pargi (nhigh)
call flush (STDOUT)
}
# Do the image average. Switch on the calculation type.
switch (calctype) {
case TY_SHORT:
call imsums (list, Memc[image], im_out, nlow, nhigh, Memc[option])
case TY_INT:
call imsumi (list, Memc[image], im_out, nlow, nhigh, Memc[option])
case TY_LONG:
call imsuml (list, Memc[image], im_out, nlow, nhigh, Memc[option])
case TY_REAL:
call imsumr (list, Memc[image], im_out, nlow, nhigh, Memc[option])
case TY_DOUBLE:
call imsumd (list, Memc[image], im_out, nlow, nhigh, Memc[option])
default:
call imsumr (list, Memc[image], im_out, nlow, nhigh, Memc[option])
}
call imunmap (im_out)
call imunmap (im_in)
# Set the header parameters.
call imtrew (list)
call imsum_hparam (list, Memc[image], Memc[hparams], Memc[option])
call imtclose (list)
call sfree (sp)
end
# IMSUM_SET -- Determine the output image pixel type and the calculation
# datatype. The default pixel types are based on the highest arithmetic
# precendence of the input images.
define NTYPES 5
procedure imsum_set (list, pixtype, calctype)
int list # List of input images
int pixtype # Pixel datatype of output image
int calctype # Pixel datatype for calculations
int i, j, nimages, max_type
pointer sp, str, im1, im2
int imtgetim(), imtlen()
bool xt_imleneq()
pointer immap()
errchk immap, imunmap
begin
call smark (sp)
call salloc (str, SZ_LINE, TY_CHAR)
# Determine maximum precedence datatype.
# Also check that the images are the same dimension and size.
nimages = imtlen (list)
j = imtgetim (list, Memc[str], SZ_LINE)
im1 = immap (Memc[str], READ_ONLY, 0)
max_type = IM_PIXTYPE (im1)
do i = 2, nimages {
j = imtgetim (list, Memc[str], SZ_LINE)
im2 = immap (Memc[str], READ_ONLY, 0)
if ((IM_NDIM(im1) != IM_NDIM(im2)) || !xt_imleneq (im1, im2)) {
call imunmap (im1)
call imunmap (im2)
call error (0, "Images have different dimensions or sizes")
}
switch (IM_PIXTYPE (im2)) {
case TY_SHORT:
if (max_type == TY_USHORT)
max_type = TY_INT
case TY_USHORT:
if (max_type == TY_SHORT)
max_type = TY_INT
case TY_INT:
if (max_type == TY_USHORT || max_type == TY_SHORT)
max_type = IM_PIXTYPE (im2)
case TY_LONG:
if (max_type == TY_USHORT || max_type == TY_SHORT ||
max_type == TY_INT)
max_type = IM_PIXTYPE (im2)
case TY_REAL:
if (max_type != TY_DOUBLE)
max_type = IM_PIXTYPE (im2)
case TY_DOUBLE:
max_type = IM_PIXTYPE (im2)
default:
}
call imunmap (im2)
}
call imunmap (im1)
call imtrew (list)
# Set calculation datatype.
call clgstr ("calctype", Memc[str], SZ_LINE)
switch (Memc[str]) {
case EOS:
calctype = max_type
case 's':
calctype = TY_SHORT
case 'i':
calctype = TY_INT
case 'l':
calctype = TY_LONG
case 'r':
calctype = TY_REAL
case 'd':
calctype = TY_DOUBLE
default:
call error (0, "Unrecognized datatype")
}
# Set output pixel datatype.
call clgstr ("pixtype", Memc[str], SZ_LINE)
switch (Memc[str]) {
case EOS:
pixtype = calctype
case 'u':
pixtype = TY_USHORT
case 's':
pixtype = TY_SHORT
case 'i':
pixtype = TY_INT
case 'l':
pixtype = TY_LONG
case 'r':
pixtype = TY_REAL
case 'd':
pixtype = TY_DOUBLE
default:
call error (0, "Unrecognized datatype")
}
call sfree (sp)
end
# IMSUM_HPARM -- Arithmetic on image header parameters.
#
# This program is limited by a lack of a rewind procedure for the image
# header fields list. Thus, a static array of field names is used
# to require only one pass through the list and the images.
define NFIELDS 10 # Maximum number of fields allowed.
procedure imsum_hparam (list, output, hparams, option)
int list # List of input images.
char output[ARB] # Output image
char hparams[ARB] # List of header parameters
char option[ARB] # Sum option
int i, nfields, flist
pointer sp, field, dvals, image, in, out
int imofnlu(), imgnfn(), imtgetim(), imtlen()
bool strne(), streq()
double imgetd()
pointer immap()
errchk immap, imofnlu, imgetd, imputd, imunmap
begin
# Return if median.
if (strne (option, "average") && strne (option, "sum"))
return
# Allocate memory.
call smark (sp)
call salloc (field, NFIELDS*SZ_FNAME, TY_CHAR)
call salloc (dvals, NFIELDS, TY_DOUBLE)
call salloc (image, SZ_FNAME, TY_CHAR)
# Map the fields.
out = immap (output, READ_WRITE, 0)
flist = imofnlu (out, hparams)
i = 0
while ((i < NFIELDS) &&
(imgnfn (flist, Memc[field+i*SZ_FNAME], SZ_FNAME) != EOF))
i = i + 1
call imcfnl (flist)
# Accumulate values from each image.
nfields = i
call aclrd (Memd[dvals], nfields)
while (imtgetim (list, Memc[image], SZ_FNAME) != EOF) {
in = immap (Memc[image], READ_ONLY, 0)
do i = 1, nfields
Memd[dvals+i-1] = Memd[dvals+i-1] +
imgetd (in, Memc[field+(i-1)*SZ_FNAME])
call imunmap (in)
}
# Output the sums or average.
if (streq (option, "average")) {
i = imtlen (list)
call adivkd (Memd[dvals], double (i), Memd[dvals], nfields)
}
do i = 1, nfields
call imputd (out, Memc[field+(i-1)*SZ_FNAME], Memd[dvals+i-1])
call imunmap (out)
call sfree (sp)
end
|
