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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <imhdr.h>
include <mwset.h>
define NTYPES 7
# T_IMSTACK -- Stack images into a single image of higher dimension.
procedure t_imstack ()
int i, j, npix, list, pdim, lmax, lindex
int axno[IM_MAXDIM], axval[IM_MAXDIM]
long line_in[IM_MAXDIM], line_out[IM_MAXDIM]
pointer sp, input, output, in, out, buf_in, buf_out, mwin, mwout
bool envgetb()
int imtopenp(), imtgetim(), imtlen()
int imgnls(), imgnli(), imgnll(), imgnlr(), imgnld(), imgnlx()
int impnls(), impnli(), impnll(), impnlr(), impnld(), impnlx()
int mw_stati()
pointer immap(), mw_open(), mw_openim()
begin
call smark (sp)
call salloc (input, SZ_FNAME, TY_CHAR)
call salloc (output, SZ_FNAME, TY_CHAR)
# Get the input images and the output image.
list = imtopenp ("images")
call clgstr ("output", Memc[output], SZ_FNAME)
# Add each input image to the output image.
i = 0
while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
i = i + 1
in = immap (Memc[input], READ_ONLY, 0)
# For the first input image map the output image as a copy
# and increment the dimension. Set the output line counter.
if (i == 1) {
out = immap (Memc[output], NEW_COPY, in)
call isk_new_image (out)
IM_NDIM(out) = IM_NDIM(out) + 1
IM_LEN(out, IM_NDIM(out)) = imtlen (list)
npix = IM_LEN(out, 1)
call amovkl (long(1), line_out, IM_MAXDIM)
}
# Check next input image for consistency with the output image.
if (IM_NDIM(in) != IM_NDIM(out) - 1)
call error (0, "Input images not consistent")
do j = 1, IM_NDIM(in) {
if (IM_LEN(in, j) != IM_LEN(out, j))
call error (0, "Input images not consistent")
}
# Copy the input lines from the image to the next lines of
# the output image. Switch on the output data type to optimize
# IMIO.
call amovkl (long(1), line_in, IM_MAXDIM)
switch (IM_PIXTYPE (out)) {
case TY_SHORT:
while (imgnls (in, buf_in, line_in) != EOF) {
if (impnls (out, buf_out, line_out) == EOF)
call error (0, "Error writing output image")
call amovs (Mems[buf_in], Mems[buf_out], npix)
}
case TY_INT:
while (imgnli (in, buf_in, line_in) != EOF) {
if (impnli (out, buf_out, line_out) == EOF)
call error (0, "Error writing output image")
call amovi (Memi[buf_in], Memi[buf_out], npix)
}
case TY_USHORT, TY_LONG:
while (imgnll (in, buf_in, line_in) != EOF) {
if (impnll (out, buf_out, line_out) == EOF)
call error (0, "Error writing output image")
call amovl (Meml[buf_in], Meml[buf_out], npix)
}
case TY_REAL:
while (imgnlr (in, buf_in, line_in) != EOF) {
if (impnlr (out, buf_out, line_out) == EOF)
call error (0, "Error writing output image")
call amovr (Memr[buf_in], Memr[buf_out], npix)
}
case TY_DOUBLE:
while (imgnld (in, buf_in, line_in) != EOF) {
if (impnld (out, buf_out, line_out) == EOF)
call error (0, "Error writing output image")
call amovd (Memd[buf_in], Memd[buf_out], npix)
}
case TY_COMPLEX:
while (imgnlx (in, buf_in, line_in) != EOF) {
if (impnlx (out, buf_out, line_out) == EOF)
call error (0, "Error writing output image")
call amovx (Memx[buf_in], Memx[buf_out], npix)
}
default:
while (imgnlr (in, buf_in, line_in) != EOF) {
if (impnlr (out, buf_out, line_out) == EOF)
call error (0, "Error writing output image")
call amovr (Memr[buf_in], Memr[buf_out], npix)
}
}
# Update the wcs. The output image will inherit the wcs of
# the first input image. The new axis will be assigned the
# identity transformation if wcsdim of the original image is
# less than the number of dimensions in the stacked image.
if ((i == 1) && (! envgetb ("nowcs"))) {
mwin = mw_openim (in)
pdim = mw_stati (mwin, MW_NPHYSDIM)
call mw_gaxmap (mwin, axno, axval, pdim)
lmax = 0
lindex = 0
do j = 1, pdim {
if (axno[j] <= lmax)
next
lmax = axno[j]
lindex = j
}
if (lindex < pdim) {
axno[pdim] = lmax + 1
axval[pdim] = 0
call mw_saxmap (mwin, axno, axval, pdim)
call mw_saveim (mwin, out)
} else {
mwout = mw_open (NULL, pdim + 1)
call isk_wcs (mwin, mwout, IM_NDIM(out))
call mw_saveim (mwout, out)
call mw_close (mwout)
}
call mw_close (mwin)
}
call imunmap (in)
}
# Finish up.
call imunmap (out)
call imtclose (list)
call sfree (sp)
end
# ISK_NEW_IMAGE -- Get a new image title and pixel type.
#
# The strings 'default' or '*' are recognized as defaulting to the original
# title or pixel datatype.
procedure isk_new_image (im)
pointer im # image descriptor
pointer sp, lbuf
int i, type_codes[NTYPES]
bool strne()
int stridx()
string types "suilrdx"
data type_codes /TY_SHORT,TY_USHORT,TY_INT,TY_LONG,TY_REAL,TY_DOUBLE,
TY_COMPLEX/
begin
call smark (sp)
call salloc (lbuf, SZ_LINE, TY_CHAR)
call clgstr ("title", Memc[lbuf], SZ_LINE)
if (strne (Memc[lbuf], "default") && strne (Memc[lbuf], "*"))
call strcpy (Memc[lbuf], IM_TITLE(im), SZ_IMTITLE)
call clgstr ("pixtype", Memc[lbuf], SZ_LINE)
if (strne (Memc[lbuf], "default") && strne (Memc[lbuf], "*")) {
i = stridx (Memc[lbuf], types)
if (i != 0)
IM_PIXTYPE(im) = type_codes[i]
}
call sfree (sp)
end
# ISK_WCS -- Update the wcs of the stacked image.
procedure isk_wcs (mwin, mwout, ndim)
pointer mwin # input wcs descriptor
pointer mwout # output wcs descriptor
int ndim # the dimension of the output image
int i, j, nin, nout, szatstr, axno[IM_MAXDIM], axval[IM_MAXDIM]
pointer sp, wcs, attribute, matin, matout, rin, rout, win, wout, atstr
int mw_stati(), itoc(), strlen()
errchk mw_newsystem()
begin
# Get the sizes of the two wcs.
nin = mw_stati (mwin, MW_NPHYSDIM)
nout = mw_stati (mwout, MW_NPHYSDIM)
szatstr = SZ_LINE
# Allocate space for the matrices and vectors.
call smark (sp)
call salloc (wcs, SZ_FNAME, TY_CHAR)
call salloc (matin, nin * nin, TY_DOUBLE)
call salloc (matout, nout * nout, TY_DOUBLE)
call salloc (rin, nin, TY_DOUBLE)
call salloc (rout, nout, TY_DOUBLE)
call salloc (win, nin, TY_DOUBLE)
call salloc (wout, nout, TY_DOUBLE)
call salloc (attribute, SZ_FNAME, TY_CHAR)
call malloc (atstr, szatstr, TY_CHAR)
# Set the system name.
call mw_gsystem (mwin, Memc[wcs], SZ_FNAME)
iferr (call mw_newsystem (mwout, Memc[wcs], nout))
call mw_ssystem (mwout, Memc[wcs])
# Set the lterm.
call mw_gltermd (mwin, Memd[matin], Memd[rin], nin)
call aclrd (Memd[rout], nout)
call amovd (Memd[rin], Memd[rout], nin)
call mw_mkidmd [Memd[matout], nout)
call isk_mcopy (Memd[matin], nin, Memd[matout], nout)
call mw_sltermd (mwout, Memd[matout], Memd[rout], nout)
# Set the wterm.
call mw_gwtermd (mwin, Memd[rin], Memd[win], Memd[matin], nin)
call aclrd (Memd[rout], nout)
call amovd (Memd[rin], Memd[rout], nin)
call aclrd (Memd[wout], nout)
call amovd (Memd[win], Memd[wout], nin)
call mw_mkidmd [Memd[matout], nout)
call isk_mcopy (Memd[matin], nin, Memd[matout], nout)
call mw_swtermd (mwout, Memd[rout], Memd[wout], Memd[matout], nout)
# Set the axis map.
call mw_gaxmap (mwin, axno, axval, nin)
do i = nin + 1, nout {
axno[i] = ndim
axval[i] = 0
}
call mw_saxmap (mwout, axno, axval, nout)
# Get the axis list and copy the old attribute list for each axis.
do i = 1, nin {
iferr (call mw_gwattrs (mwin, i, "wtype", Memc[atstr], szatstr))
call strcpy ("linear", Memc[atstr], szatstr)
call mw_swtype (mwout, i, 1, Memc[atstr], "")
for (j = 1; ; j = j + 1) {
if (itoc (j, Memc[attribute], SZ_FNAME) <= 0)
Memc[attribute] = EOS
repeat {
iferr (call mw_gwattrs (mwin, i, Memc[attribute],
Memc[atstr], szatstr))
Memc[atstr] = EOS
if (strlen (Memc[atstr]) < szatstr)
break
szatstr = szatstr + SZ_LINE
call realloc (atstr, szatstr, TY_CHAR)
}
if (Memc[atstr] == EOS)
break
call mw_swattrs (mwout, i, Memc[attribute], Memc[atstr])
}
}
# Set the default attributes for the new axes.
do i = nin + 1, nout
call mw_swtype (mwout, i, 1, "linear", "")
call mfree (atstr, TY_CHAR)
call sfree (sp)
end
# ISK_MCOPY -- Copy a smaller 2d matrix into a larger one.
procedure isk_mcopy (matin, nin, matout, nout)
double matin[nin,nin] # the input matrix
int nin # size of the input matrix
double matout[nout,nout] # the input matrix
int nout # size of the output matrix
int i,j
begin
do i = 1, nin {
do j = 1, nin
matout[j,i] = matin[j,i]
}
end
|
