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
|
# AP_GROW_REGIONS -- Perform region growing around a rejected pixel.
int procedure ap_grow_regions (skypix, coords, wgt, nskypix, sky_zero,
index, snx, sny, rgrow, sumpx, sumsqpx, sumcbpx)
real skypix[ARB] # sky pixels
int coords[ARB] # sky cordinates
real wgt[ARB] # weights
int nskypix # total number of sky pixels
real sky_zero # sky zero point for moment analysis
int index # index of pixel to be rejected
int snx, sny # size of sky subraster
real rgrow # region growing radius
double sumpx, sumsqpx, sumcbpx # sum and sum of squares of sky pixels
double dsky
int j, k, ixc, iyc, xmin, xmax, ymin, ymax, cstart, c, nreject
real rgrow2, r2, d
begin
# Find the center of the region to be rejected.
ixc = mod (coords[index], snx)
if (ixc == 0)
ixc = snx
iyc = (coords[index] - ixc) / snx + 1
# Define the region to be searched.
rgrow2 = rgrow ** 2
ymin = max (1, int (iyc - rgrow))
ymax = min (sny, int (iyc + rgrow))
xmin = max (1, int (ixc - rgrow))
xmax = min (snx, int (ixc + rgrow))
if (ymin <= iyc)
cstart = min (nskypix, max (1, index - int (rgrow) + snx *
(ymin - iyc)))
else
cstart = index
# Reject the pixels.
nreject = 0
do j = ymin, ymax {
d = rgrow2 - (j - iyc) ** 2
if (d <= 0.0)
d = 0.0
else
d = sqrt (d)
do k = max (xmin, int (ixc - d)), min (xmax, int (ixc + d)) {
c = k + (j - 1) * snx
while (coords[cstart] < c && cstart < nskypix)
cstart = cstart + 1
r2 = (k - ixc) ** 2 + (j - iyc) ** 2
if (r2 <= rgrow2 && c == coords[cstart] && wgt[cstart] > 0.0) {
dsky = skypix[cstart] - sky_zero
sumpx = sumpx - dsky
sumsqpx = sumsqpx - dsky ** 2
sumcbpx = sumcbpx - dsky ** 3
nreject = nreject + 1
wgt[cstart] = 0.0
}
}
}
return (nreject)
end
|
