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
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <fset.h>
include <imhdr.h>
include <math/curfit.h>
include <pkg/gtools.h>
include "curfit.h"
define VERBOSE_OUTPUT 1
define LIST_OUTPUT 2
define DEFAULT_OUTPUT 3
define CF_UNIFORM 1
define CF_USER 2
define CF_STATISTICAL 3
define CF_INSTRUMENTAL 4
# CF_FIT -- Called once for each curve to be fit.
$for (rd)
procedure cf_fit$t (ic, gt, x, y, wts, nvalues, nmax, device, interactive, ofmt,
power)
pointer ic # ICFIT pointer
pointer gt # Graphics tools pointer
PIXEL x[nmax] # X data values
PIXEL y[nmax] # Y data values
PIXEL wts[nmax] # Weights
int nvalues # Number of data points
int nmax # Maximum number of data points
char device[SZ_FNAME] # Output graphics device
int interactive # Fit curve interactively?
int ofmt # Type of output listing
bool power # Convert coeff to power series?
int ncoeff, i
PIXEL xmin, xmax
pointer sp, gp, cv, coeff, tty
pointer gopen(), ttyodes()
int fstati(), $tcvstati()
begin
# Determine data range and set up curve fitting limits.
call alim$t (x, nvalues, xmin, xmax)
call ic_putr (ic, "xmin", real (xmin))
call ic_putr (ic, "xmax", real (xmax))
if (interactive == YES) {
gp = gopen (device, NEW_FILE, STDGRAPH)
call icg_fit$t (ic, gp, "cursor", gt, cv, x, y, wts, nvalues)
call gclose (gp)
} else
# Do fit non-interactively
call ic_fit$t (ic, cv, x, y, wts, nvalues, YES, YES, YES, YES)
# Output answers to STDOUT
if (ofmt != LIST_OUTPUT) {
if (fstati (STDOUT, F_REDIR) == NO) {
tty = ttyodes ("terminal")
call ttyclear (STDOUT, tty)
call ttycdes (tty)
}
#call ic_show (ic, "STDOUT", gt)
call ic_vshow$t (ic, "STDOUT", cv, x, y, wts, nvalues, gt)
if (ofmt == VERBOSE_OUTPUT) {
call printf (
"\n# \t X \t Yc \t Y \t W\n")
call cf_listxy$t (cv, x, y, wts, nvalues)
}
} else
call cf_listxy$t (cv, x, y, wts, nvalues)
# Convert coefficients if requested for legendre or chebyshev
if (power && ofmt != LIST_OUTPUT) {
# Calculate and print coefficients
ncoeff = $tcvstati (cv, CVNCOEFF)
call smark (sp)
call salloc (coeff, ncoeff, TY_PIXEL)
call $tcvpower (cv, Mem$t[coeff], ncoeff)
call printf ("# Power series coefficients would be:\n")
call printf ("# \t\tcoefficient\n")
do i = 1, ncoeff {
call printf ("# \t%d \t%14.7e\n")
call pargi (i)
call parg$t (Mem$t[coeff+i-1])
}
call sfree (sp)
}
$if (datatype == r)
call cvfree (cv)
$else
call $tcvfree (cv)
$endif
#call ic_close$t (ic)
end
# CF_LISTXY -- Print answers to STDOUT as x,y pairs.
procedure cf_listxy$t (cv, xvals, yvals, wts, nvalues)
pointer cv # Pointer to curfit structure
int nvalues # Number of data values
PIXEL xvals[nvalues] # Array of x data values
PIXEL yvals[nvalues] # Array of y data values
PIXEL wts[nvalues] # Array of weights
int i
PIXEL $tcveval()
begin
do i = 1, nvalues {
call printf ("\t%14.7e \t%14.7e \t%14.7e \t%14.7e\n")
call parg$t (xvals[i])
call parg$t ($tcveval (cv, xvals[i]))
call parg$t (yvals[i])
call parg$t (wts[i])
}
end
# IM_PROJECTION -- Given an image section of arbitrary dimension, compute
# the projection along a single axis by taking the average over the other
# axes. We do not know about bad pixels.
procedure im_projection$t (im, x, y, w, npix, weighting, axis)
pointer im # Pointer to image header structure
PIXEL x[npix] # Index of projection vector
PIXEL y[npix] # Receives the projection vector
PIXEL w[npix] # Receives the weight vector
int weighting # Weighting of the individual points
int npix # Length of projection vector
int axis # The axis to be projected to (x=1)
int i, lastv
long v[IM_MAXDIM], nsum, totpix
pointer pix
PIXEL asum$t()
pointer imgnl$t()
errchk imgnl$t
begin
if (im == NULL)
call error (1, "Image projection operator called with null im")
if (axis < 1 || axis > IM_NDIM(im))
call error (2, "Attempt to take projection over nonexistent axis")
# Set the y projection vector
call aclr$t (y, npix)
call amovkl (long(1), v, IM_MAXDIM)
switch (axis) {
case 1:
# Since the image is read line by line, it is easy to compute the
# projection along the x-axis (axis 1). We merely sum all of the
# image lines.
while (imgnl$t (im, pix, v) != EOF)
call aadd$t (Mem$t[pix], y, y, npix)
default:
# Projecting along any other axis when reading the image line
# by line is a bit difficult to understand. Basically, the
# element 'axis' of the V vector (position of the line in the
# image) gives us the index into the appropriate element of
# y. When computing the projection over multiple dimensions,
# the same output element will be referenced repeatedly. All
# of the elmenents of the input line are summed and added into
# this output element.
for (lastv=v[axis]; imgnl$t (im, pix, v) != EOF; lastv=v[axis]) {
i = lastv
if (i <= npix)
y[i] = y[i] + asum$t (Mem$t[pix], IM_LEN(im,1))
}
}
# Now compute the number of pixels contributing to each element
# of the output vector. This is the number of pixels in the image
# divided by the length of the projection.
totpix = 1
do i = 1, IM_NDIM(im)
if (i == axis)
totpix = totpix * min (npix, IM_LEN(im,i))
else
totpix = totpix * IM_LEN(im,i)
nsum = totpix / min (npix, IM_LEN(im,axis))
# Compute the average by dividing by the number if pixels summed at
# each point.
call adivk$t (y, PIXEL (nsum), y, npix)
# Set the x and weight vectors
do i = 1, npix {
x[i] = i
switch (weighting) {
case CF_STATISTICAL:
if (y[i] > 0.0)
w[i] = 1.0 / y[i]
else if (y[i] < 0.0)
w[i] = abs (1.0 / y[i])
else
w[i] = 1.0
case CF_UNIFORM:
w[i] = 1.
default:
w[i] = 1.
}
}
end
$endfor
|
