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
|
# AP_GROW_HIST -- Procedure to reject pixels with region growing.
int procedure ap_grow_hist (skypix, coords, wgt, nskypix, index, snx, sny, hgm,
nbins, z1, z2, rgrow)
real skypix[ARB] # sky pixels
int coords[ARB] # sky pixel coordinates
real wgt[ARB] # array of weights
int nskypix # number of sky pixels
int index # index of pixel to be rejected
int snx, sny # size of sky subraster
real hgm[ARB] # histogram
int nbins # number of bins
real z1, z2 # value of first and last histogram bin
real rgrow # region growing radius
int j, k, ixc, iyc, ymin, ymax, xmin, xmax, nreject, cstart, c, bin
real dh, r2, rgrow2, d
begin
# Find the x and y coordinates of the pixel to be rejected.
ixc = mod (coords[index], snx)
if (ixc == 0)
ixc = snx
iyc = (coords[index] - ixc) / snx + 1
# Find the coordinate space to be used for regions growing.
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
# Perform the region growing.
dh = real (nbins - 1) / (z2 - z1)
rgrow2 = rgrow ** 2
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) {
nreject = nreject + 1
wgt[cstart] = 0.0
if (skypix[cstart] >= z1 && skypix[cstart] <= z2) {
bin = int ((skypix[cstart] - z1) * dh) + 1
hgm[bin] = hgm[bin] - 1.0
}
}
}
}
return (nreject)
end
# AP_GROW_HIST2 -- Procedure to reject pixels with region growing.
int procedure ap_grow_hist2 (skypix, coords, wgt, nskypix, sky_zero, index,
snx, sny, hgm, nbins, z1, z2, rgrow, sumpx, sumsqpx, sumcbpx)
real skypix[ARB] # sky pixels
int coords[ARB] # coordinates
real wgt[ARB] # array of weights
int nskypix # number of sky pixels
real sky_zero # the sky zero point for moment analysis
int index # index of pixel to be rejected
int snx, sny # size of sky subraster
real hgm[ARB] # histogram
int nbins # number of bins
real z1, z2 # value of first and last histogram bin
real rgrow # region growing radius
double sumpx # sum of sky values
double sumsqpx # sum of sky values squared
double sumcbpx # sum of sky values cubed
double dsky
int j, k, ixc, iyc, ymin, ymax, xmin, xmax, nreject, cstart, c, bin
real dh, r2, rgrow2, d
begin
# Find the coordinates of the region growing center.
ixc = mod (coords[index], snx)
if (ixc == 0)
ixc = snx
iyc = (coords[index] - ixc) / snx + 1
ymin = max (1, int (iyc - rgrow))
ymax = min (sny, int (iyc + rgrow))
xmin = max (1, int (ixc - rgrow))
xmax = min (snx, int (ixc + rgrow))
dh = real (nbins - 1) / (z2 - z1)
if (ymin <= iyc)
cstart = min (nskypix, max (1, index - int (rgrow) + snx *
(ymin - iyc)))
else
cstart = index
# Perform the region growing.
nreject = 0
rgrow2 = rgrow ** 2
do j = ymin, ymax {
d = sqrt (rgrow2 - (j - iyc) ** 2)
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) {
nreject = nreject + 1
wgt[cstart] = 0.0
dsky = skypix[cstart] - sky_zero
sumpx = sumpx - dsky
sumsqpx = sumsqpx - dsky ** 2
sumcbpx = sumcbpx - dsky ** 3
if (skypix[cstart] >= z1 && skypix[cstart] <= z2) {
bin = int ((skypix[cstart] - z1) * dh) + 1
hgm[bin] = hgm[bin] - 1.0
}
}
}
}
return (nreject)
end
|
