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
|
include <mach.h>
include "../lib/fitsky.h"
# AP_MEAN -- Procedure to calculate the mean of the sky pixel array.
int procedure ap_mean (skypix, coords, wgt, index, nskypix, snx, sny, losigma,
hisigma, rgrow, maxiter, sky_mean, sky_sigma, sky_skew, nsky,
nsky_reject)
real skypix[ARB] # unsorted array of skypixels
int coords[ARB] # coordinate array for region growing
real wgt[ARB] # the weight array for rejection
int index[ARB] # sorted array of indices
int nskypix # total number of sky pixels
int snx, sny # dimensions of the sky subraster
real losigma, hisigma # number of sky_sigma for rejection
real rgrow # radius of region growing
int maxiter # maximum number of cycles of rejection
real sky_mean # computed sky value
real sky_sigma # the computed sky sigma
real sky_skew # skewness of sky distribution
int nsky # the number of sky pixels used
int nsky_reject # the number of sky pixels rejected
double dsky, sumpx, sumsqpx, sumcbpx
int i, j, nreject
real sky_zero, dmin, dmax, locut, hicut
int ap_grow_regions()
real ap_asumr()
begin
# Intialize.
nsky = nskypix
nsky_reject = 0
sky_mean = INDEFR
sky_sigma = INDEFR
sky_skew = INDEFR
if (nskypix <= 0)
return (AP_NOSKYAREA)
# Compute the mean, sigma and skew.
sky_zero = ap_asumr (skypix, index, nskypix) / nskypix
call apfimoments (skypix, index, nskypix, sky_zero, sumpx, sumsqpx,
sumcbpx, sky_mean, sky_sigma, sky_skew)
call ap_ialimr (skypix, index, nskypix, dmin, dmax)
sky_mean = max (dmin, min (sky_mean, dmax))
# Decide whether to do the rejection cycle.
if (maxiter < 1 || (IS_INDEFR(losigma) && IS_INDEFR(hisigma)) ||
sky_sigma <= 0.0)
return (AP_OK)
# Reject points within k1 * sky_sigma of the median.
do i = 1, maxiter {
# Compute the new rejection limits.
if (IS_INDEFR(losigma))
locut = max (-MAX_REAL, dmin)
else if (i == 1)
locut = sky_mean - min (sky_mean - dmin, dmax - sky_mean,
losigma * sky_sigma)
else
locut = sky_mean - losigma * sky_sigma
if (IS_INDEFR(hisigma))
hicut = min (MAX_REAL, dmax)
else if (i == 1)
hicut = sky_mean + min (dmax - sky_mean, sky_mean - dmin,
hisigma * sky_sigma)
else
hicut = sky_mean + hisigma * sky_sigma
nreject = 0
do j = 1, nskypix {
if (wgt[index[j]] <= 0.0)
next
if (skypix[index[j]] < locut || skypix[index[j]] > hicut) {
if (rgrow > 0.0) {
nreject = ap_grow_regions (skypix, coords, wgt,
nskypix, sky_zero, index[j], snx, sny, rgrow,
sumpx, sumsqpx, sumcbpx)
} else {
dsky = (skypix[index[j]] - sky_zero)
sumpx = sumpx - dsky
sumsqpx = sumsqpx - dsky ** 2
sumcbpx = sumcbpx - dsky ** 3
wgt[index[j]] = 0.0
nreject = nreject + 1
}
}
}
# Test the number of rejected pixels.
if (nreject <= 0)
break
nsky_reject = nsky_reject + nreject
# Test that some pixels actually remain.
nsky = nskypix - nsky_reject
if (nsky <= 0)
break
# Recompute the mean, sigma and skew.
call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero,
sky_mean, sky_sigma, sky_skew)
if (sky_sigma <= 0.0)
break
sky_mean = max (dmin, min (sky_mean, dmax))
}
# Return an appropriate error code.
if (nsky == 0 || nsky_reject == nskypix) {
nsky = 0
nsky_reject = nskypix
sky_mean = INDEFR
sky_sigma = INDEFR
sky_skew = INDEFR
return (AP_NSKY_TOO_SMALL)
} else
return (AP_OK)
end
|
