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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <imhdr.h>
include <imio.h>
include "mwcs.h"
# IW_SETAXMAP -- If the reference image was opened with an image section,
# modify the Lterm to reflect the section transformation, and enable the
# axis map if any dimensional reduction was involved.
procedure iw_setaxmap (mw, im)
pointer mw #I pointer to MWCS descriptor
pointer im #I pointer to reference image
double v
pointer sp, ltv_1, ltv_2, ltm
int wcsdim, ndim, physax, i, j
int axno[MAX_DIM], axval[MAX_DIM]
int o_axno[MAX_DIM], o_axval[MAX_DIM]
int n_axno[MAX_DIM], n_axval[MAX_DIM]
begin
# If there is no section we don't need to do anything.
if (IM_SECTUSED(im) == NO)
return
call smark (sp)
ndim = IM_NPHYSDIM(im)
call salloc (ltv_1, ndim, TY_DOUBLE)
call salloc (ltv_2, ndim, TY_DOUBLE)
call salloc (ltm, ndim*ndim, TY_DOUBLE)
# The section transformation is px = VSTEP * lx + VOFF, specifying
# the transformation from logical to physical image coordinates.
# The IMIO axis map is given by j=VMAP[i], mapping logical axis I to
# physical axis J. Hence the physical to logical transformation in
# terms of IMIO units is given by lx = (1/VSTEP) * px + (-VOFF/VSTEP).
# Since the section transform forbids rotation the axes are independent.
call aclrd (Memd[ltv_1], ndim)
call aclrd (Memd[ltm], ndim * ndim)
do i = 1, ndim {
if (IM_VSTEP(im,i) == 0)
v = 1.0D0
else
v = 1.0D0 / IM_VSTEP(im,i)
Memd[ltm+(i-1)*ndim+i-1] = v
Memd[ltv_2+(i-1)] = -(IM_VOFF(im,i) * v)
}
# Enter the section transformation. This uses the axis map, but the
# transformation is defined in terms of the physical image matrix,
# which is defined by the old axis map before modification by the new
# image section. Hence we must do this step before editing the axis
# map below.
call mw_translated (mw, Memd[ltv_1], Memd[ltm], Memd[ltv_2], ndim)
# Compute the axis map for the active image section relative to the
# current physical image matrix.
do j = 1, ndim {
for (i=1; i <= IM_NDIM(im); i=i+1)
if (IM_VMAP(im,i) == j)
break
if (i > IM_NDIM(im)) {
axno[j] = 0
axval[j] = IM_VOFF(im,j)
} else {
axno[j] = i
axval[j] = 0
}
}
# Get the old axis map for the WCS. In the general case the WCS can
# have a dimension higher than the current image, i.e. if the current
# image was produced by extracting a section of an image of higher
# dimension. In such a case the WCS will have an axis map relating
# the physical axes of the current image back to the original physical
# system.
wcsdim = MI_NDIM(mw)
call mw_gaxmap (mw, o_axno, o_axval, wcsdim)
# Combine the old axis map and the axis map for the current image
# section. The old axis map physical->logical mapping maps WCS
# physical axes to logical axes, which are the physical axes of the
# current image. The axis map for the current image section maps the
# physical axes of the current image to the logical axes of the
# section. An axis removed in the WCS axis map is not visible in the
# image axno/axval computed above; the corresponding axis in the
# combined WCS axis map is unchanged. The remaining axes are subject
# to remapping by the mage axno/axval. This mapping may set any of
# the axes to a constant to further reduce the dimensionality of the
# logical system, however that does not concern us here, we just pass
# on the combined axno/axval vectors to mw_saxmap.
do i = 1, wcsdim {
if (o_axno[i] == 0) {
n_axno[i] = 0
n_axval[i] = o_axval[i]
} else {
physax = o_axno[i]
n_axno[i] = axno[physax]
n_axval[i] = axval[physax]
}
}
# Set the new axis map.
call mw_saxmap (mw, n_axno, n_axval, wcsdim)
call sfree (sp)
end
|
