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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
|
include <mach.h>
include <imhdr.h>
include <imset.h>
include "vt.h"
# IMFILT -- Apply a spatial averageing filter to an image by convolving the
# image with a filter kernel. Return the resulting image in a separate
# image file.
procedure imfilt (inim, outim, kernel, kxdim, kydim, el)
pointer inim, outim # input and output images
int kxdim, kydim # dimensions of convolution kernel
real kernel[kxdim, kydim] # convolution kernel
real el[LEN_ELSTRUCT] # limb ellipse structure
int nlines, linelength, startline
int linebuf, outline, i
int k, offset, x2semi
int extension, startpix, lastline
real p, n, lpix1, lpix2
pointer lines, tmpptr, outptr, inptr, sp
pointer impl2r(), imgl2r(), imfglexr()
errchk impl2r, imfglexr, imgl2r
begin
# Set up the pointer array on the stack.
call smark (sp)
call salloc (lines, kydim, TY_POINTER)
# Calculate the extension.
extension = kxdim / 2
offset = E_XCENTER[el] - E_XSEMIDIAMETER[el]
x2semi = 2 * E_XSEMIDIAMETER[el]
# Startpix is the x-coordinate of the beginning of the 1-D array
# we pass to the convolution vector routine. If wrong, return.
startpix = offset - extension
if (startpix <= 0) {
call printf ("convolution kernel too wide for this image\n")
return
}
# Get the dimensions of the image.
linelength = IM_LEN(inim, 1)
nlines = IM_LEN(inim, 2)
# Pointers to the input and the output images are passed to this
# subroutine by the user.
# Use imseti to set up the appropriate number of input buffers.
call imseti (inim, IM_NBUFS, kydim+1)
# Read in the necessary number of input image lines to initially
# fill all but one of the input line buffers.
# First, skip over all lines that are off the limb.
# The size of the output image is defined prior to the call
# to this subroutine, the output image is the same size as the
# input image.
startline = 0
Memi[lines] = NULL
# Skip over empty lines.
while (Memi[lines] == NULL) {
startline = startline + 1
Memi[lines] = imfglexr (inim, startline, el, extension)
}
# Fill (almost) the line buffer.
do linebuf = 1, kydim-2
Memi[lines+linebuf] = imfglexr (inim, linebuf+startline,
el, extension)
# Copy the first startline lines from the input image into the
# output image.
do outline = 1, startline + (kydim/2) {
# Put next line to output image, get the corresponding line from
# the input image.
inptr = imgl2r (inim, outline)
outptr = impl2r (outim, outline)
# Copy the input line into the ouput line. Strip sqib.
do i = 1, DIM_VTFD {
Memr[outptr+i-1] = Memr[inptr+i-1]/16.
}
}
# Do the convolution, output line by output line.
do outline = (kydim/2) + startline, nlines {
# Use ellipse parameters to determine where the limb
# intersections are.
p = (real(outline) - E_YCENTER[el])**2/E_YSEMIDIAMETER[el]**2
n = (1.0 - p) * E_XSEMIDIAMETER[el]**2
# The two limb points are:
lpix1 = int(-sqrt(abs(n)) + .5) + E_XCENTER[el]
lpix2 = int(sqrt(abs(n)) + .5) + E_XCENTER[el]
# Keep a copy of this input line around for filling outside
# the limb.
inptr = imgl2r (inim, outline)
# Scroll the buffer pointer array.
if (outline > ((kydim/2) + startline))
do i = 0, kydim - 2
Memi[lines+i] = Memi[lines+i+1]
# Get next line from input image, if it is off the limb then we
# are done.
tmpptr = imfglexr (inim, outline+((kydim/2)+1), el, extension)
if (tmpptr == NULL) {
lastline = outline
break
}
Memi[lines+kydim-1] = tmpptr
# Put next line to output image.
outptr = impl2r (outim, outline)
# Zero the output line.
call aclrr (Memr[outptr], DIM_VTFD)
# Here is the actual convolution, this is a do loop over the lines
# of the kernel, each call to acnvrs adds the convolution of a
# kernel line with an input line to the output line.
do k = 1, kydim
call acnvr (Memr[Memi[lines+k-1]+startpix], Memr[outptr+offset],
x2semi, kernel[1,k], kxdim)
# Fill outside the limb with orig data.
do i = 1, lpix1 {
Memr[outptr+i-1] = Memr[inptr+i-1]/16.
}
do i = lpix2, DIM_VTFD {
Memr[outptr+i-1] = Memr[inptr+i-1]/16.
}
# Roundoff adjustment.
do i = startpix, startpix+x2semi {
if (Memr[outptr+i-1] < 0.0)
Memr[outptr+i-1] = Memr[outptr+i-1] - .5
else
Memr[outptr+i-1] = Memr[outptr+i-1] + .5
}
} # End of do loop on outline.
# Clear the rest of the image.
do outline = lastline, DIM_VTFD {
# Put next line to output image, get the corresponding line from
# the input image.
inptr = imgl2r (inim, outline)
outptr = impl2r (outim, outline)
# Copy the input line into the ouput line. Strip sqib.
do i = 1, DIM_VTFD {
Memr[outptr+i-1] = Memr[inptr+i-1]/16.
}
}
call sfree (sp)
end
|
