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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <imhdr.h>
include <error.h>
include <mwset.h>
# T_IMTRANSPOSE -- Transpose images.
#
# The input and output images are given by image template lists. The
# number of output images must match the number of input images. Image
# sections are allowed in the input images and are ignored in the output
# images. If the input and output image names are the same then the transpose
# is performed to a temporary file which then replaces the input image.
procedure t_imtranspose ()
char imtlist1[SZ_LINE] # Input image list
char imtlist2[SZ_LINE] # Output image list
int len_blk # 1D length of transpose block
char image1[SZ_FNAME] # Input image name
char image2[SZ_FNAME] # Output image name
char imtemp[SZ_FNAME] # Temporary file
int list1, list2
pointer im1, im2, mw
bool envgetb()
int clgeti(), imtopen(), imtgetim(), imtlen()
pointer immap(), mw_openim()
begin
# Get input and output image template lists and the size of
# the transpose block.
call clgstr ("input", imtlist1, SZ_LINE)
call clgstr ("output", imtlist2, SZ_LINE)
len_blk = clgeti ("len_blk")
# Expand the input and output image lists.
list1 = imtopen (imtlist1)
list2 = imtopen (imtlist2)
if (imtlen (list1) != imtlen (list2)) {
call imtclose (list1)
call imtclose (list2)
call error (0, "Number of input and output images not the same")
}
# Do each set of input/output images.
while ((imtgetim (list1, image1, SZ_FNAME) != EOF) &&
(imtgetim (list2, image2, SZ_FNAME) != EOF)) {
call xt_mkimtemp (image1, image2, imtemp, SZ_FNAME)
im1 = immap (image1, READ_ONLY, 0)
im2 = immap (image2, NEW_COPY, im1)
# Do the transpose.
call imtranspose (im1, im2, len_blk)
# Update the image WCS to reflect the shift.
if (!envgetb ("nomwcs")) {
mw = mw_openim (im1)
call imtrmw (mw)
call mw_saveim (mw, im2)
call mw_close (mw)
}
# Unmap the input and output images.
call imunmap (im2)
call imunmap (im1)
call xt_delimtemp (image2, imtemp)
}
call imtclose (list1)
call imtclose (list2)
end
# IMTRANSPOSE -- Transpose an image.
#
# Divide the image into square blocks of size len_blk by len_blk.
# Transpose each block with a generic array transpose operator.
procedure imtranspose (im_in, im_out, len_blk)
pointer im_in # Input image descriptor
pointer im_out # Output image descriptor
int len_blk # 1D length of transpose block
int x1, x2, nx
int y1, y2, ny
pointer buf_in, buf_out
pointer imgs2s(), imps2s(), imgs2i(), imps2i(), imgs2l(), imps2l()
pointer imgs2r(), imps2r(), imgs2d(), imps2d(), imgs2x(), imps2x()
begin
# Output image is a copy of input image with dims transposed.
IM_LEN (im_out, 1) = IM_LEN (im_in, 2)
IM_LEN (im_out, 2) = IM_LEN (im_in, 1)
# Break the input image into blocks of at most len_blk by len_blk.
do x1 = 1, IM_LEN (im_in, 1), len_blk {
x2 = x1 + len_blk - 1
if (x2 > IM_LEN(im_in, 1))
x2 = IM_LEN(im_in, 1)
nx = x2 - x1 + 1
do y1 = 1, IM_LEN (im_in, 2), len_blk {
y2 = y1 + len_blk - 1
if (y2 > IM_LEN(im_in, 2))
y2 = IM_LEN(im_in, 2)
ny = y2 - y1 + 1
# Switch on the pixel type to optimize IMIO.
switch (IM_PIXTYPE (im_in)) {
case TY_SHORT:
buf_in = imgs2s (im_in, x1, x2, y1, y2)
buf_out = imps2s (im_out, y1, y2, x1, x2)
call imtr2s (Mems[buf_in], Mems[buf_out], nx, ny)
case TY_INT:
buf_in = imgs2i (im_in, x1, x2, y1, y2)
buf_out = imps2i (im_out, y1, y2, x1, x2)
call imtr2i (Memi[buf_in], Memi[buf_out], nx, ny)
case TY_USHORT, TY_LONG:
buf_in = imgs2l (im_in, x1, x2, y1, y2)
buf_out = imps2l (im_out, y1, y2, x1, x2)
call imtr2l (Meml[buf_in], Meml[buf_out], nx, ny)
case TY_REAL:
buf_in = imgs2r (im_in, x1, x2, y1, y2)
buf_out = imps2r (im_out, y1, y2, x1, x2)
call imtr2r (Memr[buf_in], Memr[buf_out], nx, ny)
case TY_DOUBLE:
buf_in = imgs2d (im_in, x1, x2, y1, y2)
buf_out = imps2d (im_out, y1, y2, x1, x2)
call imtr2d (Memd[buf_in], Memd[buf_out], nx, ny)
case TY_COMPLEX:
buf_in = imgs2x (im_in, x1, x2, y1, y2)
buf_out = imps2x (im_out, y1, y2, x1, x2)
call imtr2x (Memx[buf_in], Memx[buf_out], nx, ny)
default:
call error (0, "unknown pixel type")
}
}
}
end
define LTM Memd[ltr+(($2)-1)*pdim+($1)-1]
define NCD Memd[ncd+(($2)-1)*pdim+($1)-1]
define swap {temp=$1;$1=$2;$2=temp}
# IMTRMW -- Perform a transpose operation on the image WCS.
procedure imtrmw (mw)
pointer mw # pointer to the mwcs structure
int i, axes[IM_MAXDIM], axval[IM_MAXDIM]
int naxes, pdim, nelem, axmap, ax1, ax2, szatstr
pointer sp, ltr, ltm, ltv, cd, r, w, ncd, nr
pointer attribute1, attribute2, atstr1, atstr2, mwtmp
double temp
int mw_stati(), itoc(), strlen()
pointer mw_open()
errchk mw_gwattrs(), mw_newsystem()
begin
# Convert axis bitflags to the axis lists.
call mw_gaxlist (mw, 03B, axes, naxes)
if (naxes < 2)
return
# Get the dimensions of the wcs and turn off axis mapping.
pdim = mw_stati (mw, MW_NPHYSDIM)
nelem = pdim * pdim
axmap = mw_stati (mw, MW_USEAXMAP)
call mw_seti (mw, MW_USEAXMAP, NO)
szatstr = SZ_LINE
# Allocate working space.
call smark (sp)
call salloc (ltr, nelem, TY_DOUBLE)
call salloc (cd, nelem, TY_DOUBLE)
call salloc (r, pdim, TY_DOUBLE)
call salloc (w, pdim, TY_DOUBLE)
call salloc (ltm, nelem, TY_DOUBLE)
call salloc (ltv, pdim, TY_DOUBLE)
call salloc (ncd, nelem, TY_DOUBLE)
call salloc (nr, pdim, TY_DOUBLE)
call salloc (attribute1, SZ_FNAME, TY_CHAR)
call salloc (attribute2, SZ_FNAME, TY_CHAR)
# Get the wterm which corresponds to the original logical to
# world transformation.
call mw_gwtermd (mw, Memd[r], Memd[w], Memd[cd], pdim)
call mw_gltermd (mw, Memd[ltm], Memd[ltv], pdim)
call mwvmuld (Memd[ltm], Memd[r], Memd[nr], pdim)
call aaddd (Memd[nr], Memd[ltv], Memd[nr], pdim)
call mwinvertd (Memd[ltm], Memd[ltr], pdim)
call mwmmuld (Memd[cd], Memd[ltr], Memd[ncd], pdim)
# Define which physical axes the logical axes correspond to.
# and recompute the above wterm to take into account the transpose.
ax1 = axes[1]
ax2 = axes[2]
swap (NCD(ax1,ax1), NCD(ax2,ax2))
swap (NCD(ax1,ax2), NCD(ax2,ax1))
swap (Memd[w+ax1-1], Memd[w+ax2-1])
swap (Memd[nr+ax1-1], Memd[nr+ax2-1])
# Perform the transpose of the lterm.
call mw_mkidmd (Memd[ltr], pdim)
LTM(ax1,ax1) = 0.0d0
LTM(ax1,ax2) = 1.0d0
LTM(ax2,ax1) = 1.0d0
LTM(ax2,ax2) = 0.0d0
call aclrd (Memd[ltv], pdim)
call aclrd (Memd[r], pdim)
call mw_translated (mw, Memd[ltv], Memd[ltr], Memd[r], pdim)
# Get the new lterm, recompute the wterm, and store it.
call mw_gltermd (mw, Memd[ltm], Memd[ltv], pdim)
call mwmmuld (Memd[ncd], Memd[ltm], Memd[cd], pdim)
call mwinvertd (Memd[ltm], Memd[ltr], pdim)
call asubd (Memd[nr], Memd[ltv], Memd[r], pdim)
call mwvmuld (Memd[ltr], Memd[r], Memd[nr], pdim)
call mw_swtermd (mw, Memd[nr], Memd[w], Memd[cd], pdim)
# Make a new temporary wcs and set the system name.
mwtmp = mw_open (NULL, pdim)
call mw_gsystem (mw, Memc[attribute1], SZ_FNAME)
iferr (call mw_newsystem (mwtmp, Memc[attribute1], pdim))
call mw_ssystem (mwtmp, Memc[attribute1])
# Copy the wterm and the lterm to it.
call mw_gwtermd (mw, Memd[r], Memd[w], Memd[ltr], pdim)
call mw_swtermd (mwtmp, Memd[r], Memd[w], Memd[ltr], pdim)
call mw_gltermd (mw, Memd[ltr], Memd[r], pdim)
call mw_sltermd (mwtmp, Memd[ltr], Memd[r], pdim)
# Set the axis map and the axis types.
call mw_gaxmap (mw, axes, axval, pdim)
call mw_saxmap (mwtmp, axes, axval, pdim)
iferr (call mw_gwattrs (mw, ax1, "wtype", Memc[attribute1], SZ_FNAME))
call strcpy ("linear", Memc[attribute1], SZ_FNAME)
iferr (call mw_gwattrs (mw, ax2, "wtype", Memc[attribute2], SZ_FNAME))
call strcpy ("linear", Memc[attribute2], SZ_FNAME)
call mw_swtype (mwtmp, ax1, 1, Memc[attribute2], "")
call mw_swtype (mwtmp, ax2, 1, Memc[attribute1], "")
# Copy the axis attributes.
call malloc (atstr1, szatstr, TY_CHAR)
call malloc (atstr2, szatstr, TY_CHAR)
for (i = 1; ; i = i + 1) {
if (itoc (i, Memc[attribute1], SZ_FNAME) <= 0)
Memc[attribute1] = EOS
if (itoc (i, Memc[attribute2], SZ_FNAME) <= 0)
Memc[attribute2] = EOS
repeat {
iferr (call mw_gwattrs (mw, ax1, Memc[attribute1],
Memc[atstr1], szatstr))
Memc[atstr1] = EOS
iferr (call mw_gwattrs (mw, ax2, Memc[attribute2],
Memc[atstr2], szatstr))
Memc[atstr2] = EOS
if ((strlen (Memc[atstr1]) < szatstr) &&
(strlen (Memc[atstr2]) < szatstr))
break
szatstr = szatstr + SZ_LINE
call realloc (atstr1, szatstr, TY_CHAR)
call realloc (atstr2, szatstr, TY_CHAR)
}
if ((Memc[atstr1] == EOS) && (Memc[atstr2] == EOS))
break
if (Memc[atstr2] != EOS)
call mw_swattrs (mwtmp, ax1, Memc[attribute2], Memc[atstr2])
if (Memc[atstr1] != EOS)
call mw_swattrs (mwtmp, ax2, Memc[attribute1], Memc[atstr1])
}
call mfree (atstr1, TY_CHAR)
call mfree (atstr2, TY_CHAR)
call mw_close (mw)
# Delete the old wcs and set equal to the new one.
call sfree (sp)
mw = mwtmp
call mw_seti (mw, MW_USEAXMAP, axmap)
end
|
