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
|
include <error.h>
include <imhdr.h>
# T_SCOORDS -- Set sampled coordinates in spectra.
# This task is currently limited to 1D spectra.
# It reads a text file of spectral coordinates and sets the WCS.
procedure t_scoords ()
pointer speclist # List of spectrum image names
pointer coordlist # List of coordinate file names
pointer label # Coordinate axis label
pointer units # Coordinate axis units
int n, fd, open(), fscan(), nscan()
int imtopenp, imtlen(), imtgetim(), clpopnu(), fntlenb(), fntgfnb()
bool verbose, clgetb()
pointer sp, spec, coords, values, im, tmp, immap()
errchk immap, open, scoords
begin
call smark (sp)
call salloc (spec, SZ_FNAME, TY_CHAR)
call salloc (coords, SZ_FNAME, TY_CHAR)
call salloc (label, SZ_FNAME, TY_CHAR)
call salloc (units, SZ_FNAME, TY_CHAR)
# Get task parameters.
speclist = imtopenp ("images")
coordlist = clpopnu ("coords")
call clgstr ("label", Memc[label], SZ_FNAME)
call clgstr ("units", Memc[units], SZ_FNAME)
verbose = clgetb ("verbose")
# Check for match between image and coordinate lists.
if (fntlenb (coordlist) > 1 && fntlenb (coordlist) != imtlen (speclist))
call error (1, "Image and coordinate lists do not match")
# Loop through spectrum list.
while (imtgetim (speclist, Memc[spec], SZ_FNAME) != EOF) {
if (fntgfnb (coordlist, Memc[coords], SZ_FNAME) == EOF)
;
iferr {
im = NULL
fd = NULL
# Open the image.
tmp = immap (Memc[spec], READ_WRITE, 0); im = tmp
# Get the coordinate values.
tmp = open (Memc[coords], READ_ONLY, TEXT_FILE); fd = tmp
call salloc (values, IM_LEN(im,1)+1, TY_DOUBLE)
n = 0
while (fscan(fd) != EOF) {
call gargd (Memd[values+n])
if (nscan() == 1)
n = n + 1
if (n > IM_LEN(im,1))
break
}
if (n != IM_LEN(im,1))
call error (1, "Wrong number of coordinate values in file")
# Create the WCS
if (verbose) {
call printf ("SCOORDS: ")
call printf (
"Setting coordinates for %s from coordinate file %s.\n")
call pargstr (Memc[spec])
call pargstr (Memc[coords])
}
call scoords (im, Memc[label], Memc[units], Memd[values])
} then
call erract (EA_WARN)
# Close files.
if (im != NULL)
call imunmap (im)
if (fd != NULL)
call close (fd)
}
call imtclose (speclist)
call fntclsb (coordlist)
call sfree (sp)
end
# SCOORDS -- Make a multispec pixel array coordinate system.
# This is currently limited to 1D spectra.
procedure scoords (im, label, units, waves)
pointer im #I Imageio I/O pointer (must be 1D image)
char label[ARB] #I Axis label (e.g. "Wavelength")
char units[ARB] #I Axis units (e.g. "Angstroms")
double waves[ARB] #I Array of dispersion coordinates
int i, n, fd, stropen()
double dw
pointer sp, coeffs, mw, mw_open()
int axes[6]
data axes/1, 2, 1, 0, 0, 0/
errchk mw_open, mw_saveim, stropen
begin
call smark (sp)
if (IM_NDIM(im) != 1)
call error (1, "scoords: image must be one dimensional")
# Initialize the MWCS.
mw = mw_open (NULL, 2)
call mw_newsystem (mw, "multispec", 2)
call mw_swtype (mw, axes, 2, "multispec", "")
if (label[1] != EOS)
call mw_swattrs (mw, 1, "label", label)
if (units[1] != EOS)
call mw_swattrs (mw, 1, "units", label)
call mw_saxmap (mw, axes[3], axes[4], 2)
# Setup multispec coefficient string.
n = IM_LEN(im,1)
i = 20 * (n + 6)
call salloc (coeffs, i, TY_CHAR)
call aclrc (Memc[coeffs], i)
fd = stropen (Memc[coeffs], i, NEW_FILE)
# Set the common attribute parameters.
dw = (waves[n] - waves[1]) / (n - 1)
call fprintf (fd, "%d %d %d %g %g %d %g %g %g ")
call pargi (1) # Aperture number
call pargi (1) # Beam number
call pargi (2) # Dispersion type (2=non-linear)
call pargd (waves[1]) # Starting coordinate
call pargd (dw) # Average dispersion
call pargi (n) # Number of pixels
call pargd (0D0) # Redshift
call pargd (INDEFD) # Aperture limit
call pargd (INDEFD) # Aperture limit
# Set the general non-linear function parameters.
call fprintf (fd, "%g %g %d ")
call pargd (1D0) # Function weight
call pargd (0D0) # Zero point shift
call pargi (5) # Function type (5=pixel array)
# Set the pixel array function values.
call fprintf (fd, "%d")
call pargi (n) # Number of pixels
do i = 1, n {
if (i > 2) {
if ((waves[i]-waves[i-1]) * dw <= 0) {
call strclose (fd)
call mw_close (mw)
call sfree (sp)
call error (1, "Coordinates are not monotonic")
}
}
call fprintf (fd, " %g")
call pargd (waves[i]) # Coordinates
}
# Write the attribute.
call strclose (fd)
call mw_swattrs (mw, 2, "spec1", Memc[coeffs])
# Store the WCS in the image header.
call mw_saveim (mw, im)
call mw_close (mw)
call sfree (sp)
end
|
