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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <error.h>
include <syserr.h>
include <imhdr.h>
include <imio.h>
include "mwcs.h"
include "imwcs.h"
# MW_LOADIM -- Load a MWCS object saved in an image header in FITS format
# into an MWCS descriptor. Note that the MWCS descriptor is allocated
# if the input is NULL. This is to allow the WCS cards to be read to
# determine the WCS dimensionality.
procedure mw_loadim (mw, im)
pointer mw #U pointer to MWCS descriptor
pointer im #I pointer to image header
bool have_wcs
int ndim, i, j, ea_type
int axno[MAX_DIM], axval[MAX_DIM]
double maxval
pointer sp, sysname, iw, ct, wp, cp, bufp, ip
int mw_allocd(), mw_refstr(), ctoi(), envgeti()
pointer iw_rfits(), iw_findcard(), iw_gbigfits(), mw_open()
errchk iw_rfits, mw_allocd, mw_newsystem, mw_swtype, iw_enterwcs, mw_saxmap
errchk mw_open
string s_physical "physical"
define axerr_ 91
define axinit_ 92
begin
call smark (sp)
call salloc (sysname, SZ_FNAME, TY_CHAR)
# Read the FITS image header into an IMWCS descriptor.
iw = iw_rfits (mw, im, RF_REFERENCE)
if (mw == NULL) {
ndim = max (IW_NDIM(iw), IM_NPHYSDIM(im))
mw = mw_open (NULL, ndim)
}
ndim = IW_NDIM(iw)
# Initialize the MWCS descriptor from the IMWCS descriptor.
# Free any storage associated with the old descriptor.
# Start with any still allocated CTRAN descriptors.
do i = 1, MAX_CTRAN {
ct = MI_CTRAN(mw,i)
if (ct != NULL)
iferr (call mw_ctfree (ct))
call erract (EA_WARN)
}
# Free the old string and data buffers.
if (MI_SBUF(mw) != NULL)
call mfree (MI_SBUF(mw), TY_CHAR)
if (MI_DBUF(mw) != NULL)
call mfree (MI_DBUF(mw), TY_DOUBLE)
# Initialize the new descriptor.
call aclri (Memi[mw], LEN_MWCS)
MI_MAGIC(mw) = MWCS_MAGIC
MI_REFIM(mw) = im
MI_NDIM(mw) = ndim
MI_LTV(mw) = mw_allocd (mw, ndim)
MI_LTM(mw) = mw_allocd (mw, ndim * ndim)
# Set the Lterm. Set axes with no LTM scales to unit scales.
# Issue a warning by default but use "wcs_matrix_err" to allow
# setting other error actions.
call amovd (IW_LTV(iw,1), D(mw,MI_LTV(mw)), ndim)
if (iw_findcard (iw, TY_LTM, ERR, 0) != NULL) {
do i = 1, ndim {
maxval = 0D0
do j = 1, ndim {
D(mw,MI_LTM(mw)+(j-1)*ndim+(i-1)) = IW_LTM(iw,i,j)
maxval = max (maxval, abs (IW_LTM(iw,i,j)))
}
if (maxval == 0D0) {
iferr (ea_type = envgeti ("wcs_matrix_err"))
ea_type = EA_WARN
iferr {
switch (ea_type) {
case EA_FATAL, EA_ERROR:
call sprintf (Memc[sysname], SZ_FNAME,
"LTM keywords for axis %d undefined")
call pargi (i)
call error (SYS_MWMISSAX, Memc[sysname])
case EA_WARN:
IW_LTM(iw,i,i) = 1D0
D(mw,MI_LTM(mw)+(i-1)*ndim+(i-1)) = IW_LTM(iw,i,i)
call sprintf (Memc[sysname], SZ_FNAME,
"setting LTM%d_%d to %.4g")
call pargi (i)
call pargi (i)
call pargd (IW_LTM(iw,i,i))
call error (SYS_MWMISSAX, Memc[sysname])
default:
IW_LTM(iw,i,i) = 1D0
D(mw,MI_LTM(mw)+(i-1)*ndim+(i-1)) = IW_LTM(iw,i,i)
}
} then
call erract (ea_type)
}
}
} else
call mw_mkidmd (D(mw,MI_LTM(mw)), ndim)
# Set up the builtin world systems "physical" and "logical".
# Both are linear systems. The physical system is a unitary
# transformation (since world systems are defined relative to
# the physical system), and the logical system has the Lterm
# for its linear term. No wcs attributes other than wtype are
# defined.
# Create the physical system.
call mw_newsystem (mw, s_physical, ndim)
do i = 1, ndim
call mw_swtype (mw, i, 1, "linear", "")
# Create the logical system.
call mw_newsystem (mw, "logical", ndim)
do i = 1, ndim
call mw_swtype (mw, i, 1, "linear", "")
# Set W and CD for the logical system to point to the Lterm.
wp = MI_WCS(mw)
WCS_W(wp) = MI_LTV(mw)
WCS_CD(wp) = MI_LTM(mw)
# Did the image header specify a WCS?
have_wcs = false
do i = 1, IW_NCARDS(iw) {
cp = IW_CARD(iw,i)
switch (C_TYPE(cp)) {
case TY_CTYPE, TY_CRPIX, TY_CRVAL, TY_CD, TY_CDELT:
have_wcs = true
break
}
}
# Enter the saved WCS. We make up a system name for now, and patch
# it up later once the real name has been recalled along with the
# attributes.
if (have_wcs) {
call mw_newsystem (mw, "image", ndim)
call iw_enterwcs (mw, iw, ndim)
ifnoerr {
call mw_gwattrs (mw, 0, "system", Memc[sysname], SZ_FNAME)
} then
WCS_SYSTEM(MI_WCS(mw)) = mw_refstr (mw, Memc[sysname])
}
# Restore the saved WCS axis map if any.
if (iw_findcard (iw, TY_WAXMAP, ERR, 0) != NULL) {
bufp = iw_gbigfits (iw, TY_WAXMAP, ERR)
ip = bufp
do i = 1, ndim {
if (ctoi (Memc, ip, axno[i]) <= 0)
goto axerr_
if (ctoi (Memc, ip, axval[i]) <= 0) {
axerr_ call eprintf ("Image %s: cannot decode WAXMAP\n")
call pargstr (IM_NAME(IW_IM(iw)))
goto axinit_
}
}
call mfree (bufp, TY_CHAR)
call mw_saxmap (mw, axno, axval, ndim)
} else {
axinit_ do i = 1, ndim {
MI_AXNO(mw,i) = i
MI_AXVAL(mw,i) = 0
}
MI_USEAXMAP(mw) = NO
MI_NLOGDIM(mw) = ndim
}
# Apply the section transform, if the image was opened with an image
# section. This edits the axis map restored above, if any, and must
# be done after restoring the original WCS axis map.
call iw_setaxmap (mw, im)
# Set the default world system.
call mw_sdefwcs (mw)
call iw_cfits (iw)
call sfree (sp)
end
|
