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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <pkg/rg.h>
# RG_WTBIN -- Weighted average or median of data.
#
# The ranges are broken up into subranges of at most abs (nbin) points and a
# minimum of max (3, (abs(nbin)+1)/2) (though always at least one bin). The
# subranges are weighted averaged if nbin > 1 and medianed if nbin < 1.
# The output weights are the sum of the weights for each subrange.
# The output array must be large enough to contain the desired points.
# If the ranges are merged then the input and output arrays may be the same.
procedure rg_wtbin$t (rg, nbin, in, wtin, nin, out, wtout, nout)
pointer rg # Ranges
int nbin # Maximum points in average or median
PIXEL in[nin] # Input array
PIXEL wtin[nin] # Input weights
int nin # Number of input points
PIXEL out[ARB] # Output array
PIXEL wtout[ARB] # Output weights
int nout # Number of output points
int i, j, k, l, n, npts, ntemp, nsample
PIXEL asum$t(), amed$t()
errchk rg_pack$t
begin
# Check for a null set of ranges.
if (rg == NULL)
call error (0, "Range descriptor undefined")
# If the bin size is exactly one then move the selected input points
# to the output array.
if (abs (nbin) < 2) {
call rg_pack$t (rg, in, out)
call rg_pack$t (rg, wtin, wtout)
nout = RG_NPTS(rg)
return
}
# Determine the subranges and take the median or average.
npts = abs (nbin)
ntemp = 0
do i = 1, RG_NRGS(rg) {
nsample = 0
if (RG_X1(rg, i) > RG_X2(rg, i)) {
j = min (nin, RG_X1(rg, i))
k = max (1, RG_X2(rg, i))
while (j >= k) {
n = max (0, min (npts, j - k + 1))
if (nsample > 0 && n < max (min (npts, 3), (npts+1)/2))
break
k = k - n
nsample = nsample + 1
ntemp = ntemp + 1
wtout[ntemp] = asum$t (wtin[k + 1], n)
if (nbin > 0) {
if (wtout[ntemp] != 0.) {
out[ntemp] = 0.
do l = k + 1, k + n
out[ntemp] = out[ntemp] + in[l] * wtin[l]
out[ntemp] = out[ntemp] / wtout[ntemp]
} else {
out[ntemp] = 0.
do l = k + 1, k + n
out[ntemp] = out[ntemp] + in[l]
out[ntemp] = out[ntemp] / n
}
} else {
out[ntemp] = amed$t (in[k+1], n)
}
}
} else {
j = max (1, RG_X1(rg, i))
k = min (nin, RG_X2(rg, i))
while (j <= k) {
n = max (0, min (npts, k - j + 1))
if (nsample > 0 && n < max (min (npts, 3), (npts+1)/2))
break
nsample = nsample + 1
ntemp = ntemp + 1
wtout[ntemp] = asum$t (wtin[j], n)
if (nbin > 0) {
if (wtout[ntemp] != 0.) {
out[ntemp] = 0.
do l = j, j + n - 1
out[ntemp] = out[ntemp] + in[l] * wtin[l]
out[ntemp] = out[ntemp] / wtout[ntemp]
} else {
out[ntemp] = 0.
do l = j, j + n - 1
out[ntemp] = out[ntemp] + in[l]
out[ntemp] = out[ntemp] / n
}
} else {
out[ntemp] = amed$t (in[j], n)
}
j = j + n
}
}
}
nout = ntemp
end
|
