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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <imhdr.h>
include <error.h>
define AVG 1 # operation = arithmetic average
define SUM 2 # operation = arithmetic sum
# change to (lrdx) in future
$for (lrd)
# BLKAVG -- Block average or sum on n-dimensional images.
# For example, block averaging by a factor of 2 means that pixels 1 and 2
# are averaged to produce the first output pixel, 3 and 4 are averaged to
# produce the second output pixel, and so on. Remainder pixels at the end
# of a line, col, etc. are averaged correctly.
procedure blkav$t (im1, im2, blkfac, option)
pointer im1 # input image
pointer im2 # output image
int blkfac[IM_MAXDIM] # blocking factors
int option # block operation (average, sum, ...)
int num_oblks[IM_MAXDIM], i, count, ndim, dim, nlines_in_sum, nfull_blkx
int junk, num_ilines, num_olines, oline
long blkin_s[IM_MAXDIM], blkin_e[IM_MAXDIM]
long vin_s[IM_MAXDIM], vin_e[IM_MAXDIM], vout[IM_MAXDIM]
$if (datatype != l)
PIXEL sum
$else
real sum
$endif
pointer sp, accum_ptr, iline_ptr, oline_ptr
int blkcomp(), imggs$t(), impnl$t()
errchk imggs$t(), impnl$t()
begin
call smark (sp)
call salloc (accum_ptr, IM_LEN(im1, 1), TY_PIXEL)
# Initialize; input and output vectors, block counters.
ndim = IM_NDIM(im1)
nfull_blkx = IM_LEN(im1, 1) / blkfac[1]
blkin_s[1] = long(1)
blkin_e[1] = long(IM_LEN(im1, 1))
vin_s[1] = blkin_s[1]
vin_e[1] = blkin_e[1]
do i = 1, ndim {
num_oblks[i] = (IM_LEN(im1,i) + blkfac[i] - 1) / blkfac[i]
IM_LEN(im2, i) = num_oblks[i]
}
num_olines = 1
do i = 2, ndim
num_olines = num_olines * num_oblks[i]
call amovkl (long(1), vout, IM_MAXDIM)
# For each sequential output-image line, ...
do oline = 1, num_olines {
call aclr$t (Mem$t[accum_ptr], IM_LEN(im1, 1))
nlines_in_sum = 0
# Compute block vector; initialize dim>1 input vector.
num_ilines = blkcomp (im1, blkfac, vout, blkin_s, blkin_e,
vin_s, vin_e)
# Output pointer; note impnl$t returns vout for NEXT output line.
junk = impnl$t (im2, oline_ptr, vout)
# For all input lines mapping to current output line, ...
do i = 1, num_ilines {
# Get line from input image.
iline_ptr = imggs$t (im1, vin_s, vin_e, ndim)
# Increment general section input vector between block bounds.
do dim = 2, ndim {
if (vin_s[dim] < blkin_e[dim]) {
vin_s[dim] = vin_s[dim] + 1
vin_e[dim] = vin_s[dim]
break
} else {
vin_s[dim] = blkin_s[dim]
vin_e[dim] = vin_s[dim]
}
}
# Accumulate line into block sum. Keep track of no. of
# lines in sum so that we can compute block average later.
call aadd$t (Mem$t[iline_ptr], Mem$t[accum_ptr],
Mem$t[accum_ptr], IM_LEN(im1,1))
nlines_in_sum = nlines_in_sum + 1
}
# We now have a templine of sums; average/sum into output buffer
# first the full blocks using a vop.
if (option == AVG)
call abav$t (Mem$t[accum_ptr], Mem$t[oline_ptr], nfull_blkx,
blkfac[1])
else
call absu$t (Mem$t[accum_ptr], Mem$t[oline_ptr], nfull_blkx,
blkfac[1])
# Now average/sum the final partial block in x, if any.
if (nfull_blkx < num_oblks[1]) {
sum = 0$f
count = 0
do i = nfull_blkx * blkfac[1] + 1, IM_LEN(im1,1) {
sum = sum + Mem$t[accum_ptr+i-1]
count = count + 1
}
if (option == AVG)
Mem$t[oline_ptr+num_oblks[1]-1] = sum / count
else
Mem$t[oline_ptr+num_oblks[1]-1] = sum
}
# Block average into output line from the sum of all lines block
# averaged in X.
if (option == AVG)
call adivk$t (Mem$t[oline_ptr], PIXEL(nlines_in_sum),
Mem$t[oline_ptr], num_oblks[1])
}
call sfree (sp)
end
$endfor
|
