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
|
include <math/curfit.h>
# AP_RMMEASURE -- Procedure to compute the sums and areas inside the apertures.
procedure ap_rmmeasure (pixels, nx, ny, wx, wy, aperts, sums, areas, naperts)
real pixels[nx,ARB] # subraster pixel values
int nx, ny # dimensions of the subraster
real wx, wy # center of subraster
real aperts[ARB] # array of apertures
double sums[ARB] # array of sums
double areas[ARB] # aperture areas
int naperts # number of apertures
int i, j, k
double fctn
real apmaxsq, dy2, r2, r
begin
# Initialize.
apmaxsq = (aperts[naperts] + 0.5) ** 2
call aclrd (sums, naperts)
call aclrd (areas, naperts)
# Loop over the pixels.
do j = 1, ny {
dy2 = (j - wy) ** 2
do i = 1, nx {
r2 = (i - wx) ** 2 + dy2
if (r2 > apmaxsq)
next
r = sqrt (r2) - 0.5
do k = 1, naperts {
if (r > aperts[k])
next
fctn = max (0.0, min (1.0, aperts[k] - r))
sums[k] = sums[k] + fctn * pixels[i,j]
areas[k] = areas[k] + fctn
}
}
}
end
# AP_BRMMEASURE -- Procedure to compute the sums and areas inside the apertures.
procedure ap_brmmeasure (pixels, nx, ny, wx, wy, datamin, datamax, aperts,
sums, areas, naperts, minapert)
real pixels[nx,ARB] # subraster pixel values
int nx, ny # dimensions of the subraster
real wx, wy # center of subraster
real datamin # minimum good data value
real datamax # maximum good data value
real aperts[ARB] # array of apertures
double sums[ARB] # array of sums
double areas[ARB] # aperture areas
int naperts # number of apertures
int minapert # minimum number of apertures
int i, j, k, kindex
double fctn
real apmaxsq, dy2, r2, r, pixval
begin
# Initialize.
apmaxsq = (aperts[naperts] + 0.5) ** 2
call aclrd (sums, naperts)
call aclrd (areas, naperts)
minapert = naperts + 1
# Loop over the pixels.
do j = 1, ny {
dy2 = (j - wy) ** 2
do i = 1, nx {
r2 = (i - wx) ** 2 + dy2
if (r2 > apmaxsq)
next
r = sqrt (r2) - 0.5
kindex = naperts + 1
pixval = pixels[i,j]
do k = 1, naperts {
if (r > aperts[k])
next
kindex = min (k, kindex)
fctn = max (0.0, min (1.0, aperts[k] - r))
sums[k] = sums[k] + fctn * pixval
areas[k] = areas[k] + fctn
}
if (kindex < minapert) {
if (pixval < datamin || pixval > datamax)
minapert = kindex
}
}
}
end
|
