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
|
include <error.h>
include <imhdr.h>
include <mwset.h>
# T_WCSRESET -- Initialize the image wcs. The user can initialize the
# pre-defined "physical" or "world" coodinate systems, or a named
# user world coordinate system, for example the "multipsec" world
# coordinate system. If the image does not have a previously defined wcs
# then wcsreset will create the identify wcs.
procedure t_wcsreset ()
bool verbose
int ndim
pointer sp, imnamelist, image, wcs, system
pointer r, w, cd, ncd, nr, ltv, iltm, ltm
pointer imlist, im, mwim, mw
bool clgetb(), streq()
int imtgetim(), mw_stati()
pointer imtopen(), immap(), mw_openim(), mw_open()
errchk mw_openim()
begin
# Allocate working space.
call smark (sp)
call salloc (imnamelist, SZ_FNAME, TY_CHAR)
call salloc (image, SZ_FNAME, TY_CHAR)
call salloc (wcs, SZ_FNAME, TY_CHAR)
call salloc (system, SZ_FNAME, TY_CHAR)
# Get the parameters.
call clgstr ("image", Memc[imnamelist], SZ_FNAME)
call clgstr ("wcs", Memc[wcs], SZ_FNAME)
verbose = clgetb ("verbose")
# Loop through the list of images.
imlist = imtopen (Memc[imnamelist])
while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) {
# Remove any image section.
call imgimage (Memc[image], Memc[image], SZ_FNAME)
# Open the image.
im = immap (Memc[image], READ_WRITE, 0)
iferr {
if (verbose) {
call printf ("Initializing wcs %s for image %s\n")
call pargstr (Memc[wcs])
call pargstr (Memc[image])
}
mwim = mw_openim (im)
} then {
mwim = NULL
} else {
call mw_gsystem (mwim, Memc[system], SZ_FNAME)
}
# Reset the lterm only if the wcs is "physical".
if (streq (Memc[wcs], "physical") && mwim != NULL) {
# Allocate space for the transforms.
ndim = mw_stati (mwim, MW_NPHYSDIM)
call malloc (r, ndim * ndim, TY_DOUBLE)
call malloc (w, ndim * ndim, TY_DOUBLE)
call malloc (cd, ndim * ndim, TY_DOUBLE)
call malloc (ltm, ndim * ndim, TY_DOUBLE)
call malloc (ltv, ndim, TY_DOUBLE)
call malloc (iltm, ndim * ndim, TY_DOUBLE)
call malloc (nr, ndim * ndim, TY_DOUBLE)
call malloc (ncd, ndim * ndim, TY_DOUBLE)
call mw_gwtermd (mwim, Memd[r], Memd[w], Memd[cd], ndim)
call mw_gltermd (mwim, Memd[ltm], Memd[ltv], ndim)
call mwvmuld (Memd[ltm], Memd[r], Memd[nr], ndim)
call aaddd (Memd[nr], Memd[ltv], Memd[nr], ndim)
call mwinvertd (Memd[ltm], Memd[iltm], ndim)
call mwmmuld (Memd[cd], Memd[iltm], Memd[ncd], ndim)
call mw_swtermd (mwim, Memd[nr], Memd[w], Memd[ncd], ndim)
call wcs_terminit (Memd[ltm], Memd[ltv], ndim)
call mw_sltermd (mwim, Memd[ltm], Memd[ltv], ndim)
call mw_saveim (mwim, im)
# Free the space.
call mfree (r, TY_DOUBLE)
call mfree (w, TY_DOUBLE)
call mfree (cd, TY_DOUBLE)
call mfree (ncd, TY_DOUBLE)
call mfree (nr, TY_DOUBLE)
call mfree (ltm, TY_DOUBLE)
call mfree (ltv, TY_DOUBLE)
call mfree (iltm, TY_DOUBLE)
# Cannot replace physical system for unknown world system.
} else if (streq (Memc[wcs], "physical") && mwim == NULL) {
if (verbose) {
call printf ("\tCannot initialize wcs %s for image %s\n")
call pargstr (Memc[wcs])
call pargstr (Memc[image])
}
} else if (streq (Memc[wcs], "world") || streq (Memc[wcs],
Memc[system])) {
ndim = IM_NDIM(im)
mw = mw_open (NULL, ndim)
call mw_saveim (mw, im)
call mw_close (mw)
# The named wcs is not present.
} else {
call eprintf ("\tCannot find wcs %s\n")
call pargstr (Memc[wcs])
}
if (mwim != NULL)
call mw_close (mwim)
call imunmap (im)
}
call imtclose (imlist)
call sfree (sp)
end
# WCS_TERMINIT -- Initialize the shift term and rotation matrix.
procedure wcs_terminit (ltm, ltv, ndim)
double ltm[ndim,ndim] # the rotation matrix
double ltv[ndim] # the shift vector
int ndim # the number of dimensions
int i
begin
call aclrd (ltm, ndim * ndim)
do i = 1, ndim
ltm[i,i] = 1.0d0
call aclrd (ltv, ndim)
end
|
