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
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
include <mach.h>
include <imhdr.h>
# T_IMSCALE -- Scale an image.
#
# Compute the image mean between the upper and lower limits and
# scale the image to a new mean. The output image is of pixel type real.
procedure t_imscale ()
char input[SZ_FNAME] # Input image
char output[SZ_FNAME] # Output image
real mean # Output mean
real lower # Lower limit for mean
real upper # Upper limit for mean
bool verbose # Verbose output?
int i
long line_in[IM_MAXDIM], line_out[IM_MAXDIM]
real mean_in, scale
pointer in, out, data_in, data_out
int imgnlr(), impnlr()
real clgetr(), image_mean()
bool clgetb()
pointer immap()
begin
# Access images and set parameters.
call clgstr ("input", input, SZ_FNAME)
in = immap (input, READ_WRITE, 0)
call clgstr ("output", output, SZ_FNAME)
out = immap (output, NEW_COPY, in)
mean = clgetr ("mean")
lower = clgetr ("lower")
if (IS_INDEFR (lower))
lower = -MAX_REAL
upper = clgetr ("upper")
if (IS_INDEFR (upper))
upper = MAX_REAL
verbose = clgetb ("verbose")
# Set output pixel type to TY_REAL.
IM_PIXTYPE(out) = TY_REAL
# Find the image mean and rescaling.
mean_in = image_mean (in, lower, upper)
scale = mean / mean_in
# Create the output image.
call amovkl (long(1), line_in, IM_MAXDIM)
call amovkl (long(1), line_out, IM_MAXDIM)
# Loop through the image lines and rescale.
while (impnlr (out, data_out, line_out) != EOF) {
i = imgnlr (in, data_in, line_in)
call amulkr (Memr[data_in], scale, Memr[data_out], IM_LEN(in, 1))
}
if (verbose) {
call printf ("Task imscale:\n")
call printf (" Lower = %g\n")
call pargr (lower)
call printf (" Upper = %g\n")
call pargr (upper)
call printf (" %s: Mean = %g\n")
call pargstr (input)
call pargr (mean_in)
call printf (" Scale = %g\n")
call pargr (scale)
call printf (" %s: Mean = %g\n")
call pargstr (output)
call pargr (mean)
}
# Finish up
call imunmap (in)
call imunmap (out)
end
# IMAGE_MEAN -- Determine the mean value of an image between lower and upper.
#
# The algorithm here is a straight image average. In future this
# should be optimized with subsampling.
real procedure image_mean (im, lower, upper)
pointer im # IMIO descriptor
real lower # Low cutoff
real upper # High cutoff
int i, npix
long line[IM_MAXDIM]
real sum
pointer data, data_end
int imgnls(), imgnli(), imgnll(), imgnlr()
begin
sum = 0.
npix = 0
call amovkl (long(1), line, IM_MAXDIM)
# Loop through the image lines to compute the mean.
# Optimize IMIO for the image datatype.
switch (IM_PIXTYPE (im)) {
case TY_SHORT:
while (imgnls (im, data, line) != EOF) {
data_end = data + IM_LEN(im, 1) - 1
do i = data, data_end {
if ((Mems[i] < lower) || (Mems[i] > upper))
next
sum = sum + Mems[i]
npix = npix + 1
}
}
case TY_INT:
while (imgnli (im, data, line) != EOF) {
data_end = data + IM_LEN(im, 1) - 1
do i = data, data_end {
if ((Memi[i] < lower) || (Memi[i] > upper))
next
sum = sum + Memi[i]
npix = npix + 1
}
}
case TY_LONG:
while (imgnll (im, data, line) != EOF) {
data_end = data + IM_LEN(im, 1) - 1
do i = data, data_end {
if ((Meml[i] < lower) || (Meml[i] > upper))
next
sum = sum + Meml[i]
npix = npix + 1
}
}
default:
while (imgnlr (im, data, line) != EOF) {
data_end = data + IM_LEN(im, 1) - 1
do i = data, data_end {
if ((Memr[i] < lower) || (Memr[i] > upper))
next
sum = sum + Memr[i]
npix = npix + 1
}
}
}
return (sum / npix)
end
|
