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
|
include <imhdr.h>
include <mwset.h>
include "imtab.h"
# This file contains three routines: itb_wcs_init and the single and
# double precision routines itb_ctranr & itb_ctrand.
#
# Phil Hodge, 30-Sep-1993 Subroutines created.
# itb_wcs_init -- open wcs, etc
# This routine gets the wcs, turns axis mapping off, and initializes
# the transformation for physical or world coordinates. The dimension
# of the original image and the mapping from logical to physical axis
# numbers are returned for use by itb_ctrand and itb_ctranr.
procedure itb_wcs_init (im, wcs_type, mw, ct, ax, wcsdim)
pointer im # i: imhdr pointer for image
int wcs_type # i: wcs type
pointer mw, ct # o: mwcs pointers
int ax[IM_MAXDIM] # o: ax[i] is physical axis for logical axis i
int wcsdim # o: dimension of physical image coord system
#--
int axno[IM_MAXDIM] # axno[j] is logical axis for physical axis j
int axval[IM_MAXDIM] # axval[j] is value if axno[j] is zero
int ndim # number of "logical" axes
int i, j
pointer mw_openim(), mw_sctran()
int mw_stati()
begin
if (wcs_type == IMTAB_NO_WCS || wcs_type == IMTAB_LOGICAL) {
mw = NULL
ct = NULL
wcsdim = IM_NDIM(im)
return
}
# Get the wcs.
mw = mw_openim (im)
# Set up the transformation.
call mw_seti (mw, MW_USEAXMAP, NO)
if (wcs_type == IMTAB_PHYSICAL)
ct = mw_sctran (mw, "logical", "physical", 0)
else if (wcs_type == IMTAB_WORLD)
ct = mw_sctran (mw, "logical", "world", 0)
wcsdim = mw_stati (mw, MW_NPHYSDIM)
ndim = IM_NDIM(im)
# Get the logical axis number corresponding to each physical axis.
call mw_gaxmap (mw, axno, axval, wcsdim)
# Invert axno: get the physical axis number for each logical axis.
do i = 1, ndim # initialize
ax[i] = 0
do j = 1, wcsdim {
do i = 1, ndim {
if (axno[j] == i) {
ax[i] = j
break
}
}
}
# Check to be sure each axis was found.
do i = 1, ndim {
if (ax[i] < 1)
call error (1, "itb_mwcs_init: an axis was not found")
}
end
# itb_ctran -- translate coordinates with axis mapping = NO
# This routine translates "logical" coordinates to "physical" or "world".
# Axis mapping must have been turned off, and the mapping from logical
# to physical axes is given by the array AX: if I is a logical axis
# number, AX[I] is the corresponding physical axis number. Each element
# of the array IN must have been initialized to one by the calling routine.
# Separate single and double precision versions are included.
procedure itb_ctrand (im, ct, ax, in, incoords, outcoords, wcsdim)
pointer im # i: imhdr pointer
pointer ct # i: coordinate transformation pointer
int ax[wcsdim] # i: "logical" to "physical" mapping
double in[IM_MAXDIM] # io: copy of incoords but includes axis mapping
double incoords[wcsdim] # i: input "logical" coordinates
double outcoords[wcsdim] # o: output coordinates
int wcsdim # i: length of incoords & outcoords arrays
#--
int i
begin
if (ct == NULL) {
call amovd (incoords, outcoords, wcsdim)
return
}
# Take account of axis mapping; i is the logical axis number.
do i = 1, IM_NDIM(im)
in[ax[i]] = incoords[i]
call mw_ctrand (ct, in, outcoords, wcsdim)
end
procedure itb_ctranr (im, ct, ax, in, incoords, outcoords, wcsdim)
pointer im # i: imhdr pointer
pointer ct # i: coordinate transformation pointer
int ax[wcsdim] # i: "logical" to "physical" mapping
real in[IM_MAXDIM] # io: copy of incoords but includes axis mapping
real incoords[wcsdim] # i: input "logical" coordinates
real outcoords[wcsdim] # o: output coordinates
int wcsdim # i: length of incoords & outcoords arrays
#--
int i
begin
if (ct == NULL) {
call amovr (incoords, outcoords, wcsdim)
return
}
# Take account of axis mapping; i is the logical axis number.
do i = 1, IM_NDIM(im)
in[ax[i]] = incoords[i]
call mw_ctranr (ct, in, outcoords, wcsdim)
end
|
