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
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <syserr.h>
include <pmset.h>
include "../qpio.h"
define RLI_NEXTLINE 9998
define RLI_INITIALIZE 9999
define SZ_CODE 7
# QPIO_GETEVENTS -- Return a sequence of events sharing the same mask value
# which satisfy the current event attribute filter. The returned events will
# be only those in a rectangular subregion of the image (specified by a prior
# call to qpio_setrange) which are also visible through the current mask.
# Sequences of events are returned in storage order until the region is
# exhausted, at which time EOF is returned.
#
# NOTE - If debug statements (printfs) are placed in this code they will cause
# i/o problems at runtime due to reentrancy, since this routine is called in
# a low level FIO pseudodevice driver (QPF). This is also true of any of the
# routines called by this procedure, and of the related routine QPIO_READPIX.
int procedure qpio_gvtevents (io, o_ev, maskval, maxev, o_nev)
pointer io #I QPIO descriptor
pointer o_ev[maxev] #O receives the event struct pointers
int maskval #O receives the mask value of the events
int maxev #I max events out
int o_nev #O same as function value (nev_out|EOF)
int status
char code[SZ_CODE]
int qpx_gvs(), qpx_gvi(), qpx_gvl(), qpx_gvr(), qpx_gvd()
errchk syserrs
define err_ 91
begin
# The generic routines currently require that X,Y be the same type.
# It wouldn't be hard to remove this restriction if necessary, but
# it simplifies things and I doubt if a mixed types feature would
# be used very often.
if (IO_EVXTYPE(io) != IO_EVYTYPE(io))
goto err_
# Get the events.
switch (IO_EVXTYPE(io)) {
case TY_SHORT:
status = qpx_gvs (io, o_ev, maskval, maxev, o_nev)
case TY_INT:
status = qpx_gvi (io, o_ev, maskval, maxev, o_nev)
case TY_LONG:
status = qpx_gvl (io, o_ev, maskval, maxev, o_nev)
case TY_REAL:
status = qpx_gvr (io, o_ev, maskval, maxev, o_nev)
case TY_DOUBLE:
status = qpx_gvd (io, o_ev, maskval, maxev, o_nev)
default:
err_ call sprintf (code, SZ_CODE, "%d")
call pargi (IO_EVXTYPE(io))
call syserrs (SYS_QPINVEVT, code)
}
return (status)
end
$for (silrd)
# QPX_GV -- Internal generic code for qpio_getevents. There is one copy
# of this routine for each event coordinate datatype. The optimization
# strategy used here assumes that executing qpio_gv is much more expensive
# than building the call in qpio_getevents. This will normally be the case
# for a large event list or a complex expression, otherwise the operation
# is likely to be fast enough that it doesn't matter anyway.
int procedure qpx_gv$t (io, o_ev, maskval, maxev, o_nev)
pointer io #I QPIO descriptor
pointer o_ev[maxev] #O receives the event struct pointers
int maskval #O receives the mask value of the events
int maxev #I max events out
int o_nev #O same as function value (nev_out|EOF)
int x1, x2, y1, y2, xs, xe, ys, ye, x, y
pointer pl, rl, rp, bp, ex, ev, ev_p, bbmask, bb_bufp
bool useindex, lineio, bbused, rmused, nodata
int bb_xsize, bb_ysize, bb_xblock, bb_yblock, ii, jj
int v[NDIM], szs_event, mval, nev, evidx, evtop, temp, i
int ev_xoff, ev_yoff
pointer plr_open()
bool pl_linenotempty(), pl_sectnotempty()
int qpio_rbucket(), qpex_evaluate(), btoi(), plr_getpix()
define swap {temp=$1;$1=$2;$2=temp}
define putevent_ 91
define again_ 92
define done_ 93
define exit_ 94
begin
pl = IO_PL(io) # pixel list (region mask) descriptor
rl = IO_RL(io) # range list buffer
bp = IO_BP(io) # bucket buffer (type short)
ex = IO_EX(io) # QPEX (EAF) descriptor
# The following is executed when the first i/o is performed on a new
# region, to select the most efficient type of i/o to be performed,
# and initialize the i/o parameters for that case. The type of i/o
# to be performed depends upon whether or not an index can be used,
# and whether or not there is a region mask (RM) or bounding box (BB).
# The presence or absence of an event attribute filter (EAF) is not
# separated out as a special case, as it is quick and easy to test
# for the presence of an EAF and apply one it if it exists.
if (IO_ACTIVE(io) == NO) {
# Check for an index. We have an index if the event list is
# indexed, and the index is defined on the Y-coordinate we will
# be using for extraction.
useindex = (IO_INDEXLEN(io) == IO_NLINES(io) &&
IO_EVYOFF(io) == IO_IXYOFF(io) &&
IO_NOINDEX(io) == NO)
# Initialize the V and VN vectors.
do i = 1, NDIM {
IO_VN(io,i) = IO_VE(io,i) - IO_VS(io,i) + 1
if (IO_VN(io,i) < 0) {
swap (IO_VS(io,i), IO_VE(io,i))
IO_VN(io,i) = -IO_VN(io,i)
}
}
call amovi (IO_VS(io,1), IO_V(io,1), NDIM)
# Determine if full lines are to be accessed, and if a bounding
# box (subraster of the image) is defined.
lineio = (IO_VS(io,1) == 1 && IO_VE(io,1) == IO_NCOLS(io))
bbused = (!lineio || IO_VS(io,2) > 1 || IO_VE(io,2) < IO_NLINES(io))
# Determine if region mask data is to be used and if there is any
# data to be read.
nodata = (IO_NEVENTS(io) <= 0)
rmused = false
if (pl != NULL)
if (pl_sectnotempty (pl, IO_VS(io,1), IO_VE(io,1), NDIM))
rmused = true
else
nodata = true
# Select the optimal type of i/o to be used for extraction.
if (nodata) {
IO_IOTYPE(io) = NoDATA_NoAREA
useindex = false
bbused = false
} else if (bbused || rmused) {
if (useindex)
IO_IOTYPE(io) = INDEX_RMorBB
else
IO_IOTYPE(io) = NoINDEX_RMorBB
} else {
# If we are reading the entire image (no bounding box) and
# we are not using a mask, then there is no point in using
# indexed i/o.
IO_IOTYPE(io) = NoINDEX_NoRMorBB
useindex = false
}
# Initialize the range list data if it will be used.
if (useindex) {
# Dummy range specifying full line segment.
RLI_LEN(rl) = RL_FIRST
RLI_AXLEN(rl) = IO_NCOLS(io)
rp = rl + ((RL_FIRST - 1) * RL_LENELEM)
Memi[rp+RL_XOFF] = IO_VS(io,1)
Memi[rp+RL_NOFF] = IO_VN(io,1)
Memi[rp+RL_VOFF] = 1
IO_RLI(io) = RLI_INITIALIZE
}
# Open the mask for random access if i/o is not indexed and
# a region mask is used.
bbmask = IO_BBMASK(io)
if (bbmask != NULL)
call plr_close (bbmask)
if (IO_IOTYPE(io) == NoINDEX_RMorBB && rmused) {
bbmask = plr_open (pl, v, 0) # (v is never referenced)
call plr_setrect (bbmask, IO_VS(io,1),IO_VS(io,2),
IO_VE(io,1),IO_VE(io,2))
call plr_getlut (bbmask,
bb_bufp, bb_xsize, bb_ysize, bb_xblock, bb_yblock)
}
# Update the QPIO descriptor.
IO_LINEIO(io) = btoi(lineio)
IO_RMUSED(io) = btoi(rmused)
IO_BBUSED(io) = btoi(bbused)
IO_BBMASK(io) = bbmask
IO_EVI(io) = 1
IO_BKNO(io) = 0
IO_BKLASTEV(io) = 0
IO_ACTIVE(io) = YES
}
# Initialize event extraction parameters.
szs_event = IO_EVENTLEN(io)
maskval = 0
nev = 0
ev_xoff = IO_EVXOFF(io)
ev_yoff = IO_EVYOFF(io)
# Extract events using the most efficient type of i/o for the given
# selection critera (index, mask, BB, EAF, etc.).
again_
switch (IO_IOTYPE(io)) {
case NoDATA_NoAREA:
# We know in advance that there are no events to be returned,
# either because there is no data, or the area of the region
# mask within the bounding box is empty.
goto exit_
case NoINDEX_NoRMorBB:
# This is the simplest case; no index, region mask, or bounding
# box. Read and output all events in sequence.
# Refill the event bucket?
if (IO_EVI(io) > IO_BKLASTEV(io))
if (qpio_rbucket (io, IO_EVI(io)) == EOF)
goto exit_
# Copy out the event pointers.
ev = bp + (IO_EVI(io) - IO_BKFIRSTEV(io)) * szs_event
nev = min (maxev, IO_BKLASTEV(io) - IO_EVI(io) + 1)
do i = 1, nev {
o_ev[i] = ev
ev = ev + szs_event
}
IO_EVI(io) = IO_EVI(io) + nev
maskval = 1
case NoINDEX_RMorBB:
# Fully general selection, including any combination of bounding
# box, region mask, or EAF, but no index, either because there is
# no index for this event list, or the index is for a different Y
# attribute than the one being used for extraction.
bbused = (IO_BBUSED(io) == YES)
x1 = IO_VS(io,1); x2 = IO_VE(io,1)
y1 = IO_VS(io,2); y2 = IO_VE(io,2)
# Refill the event bucket?
while (IO_EVI(io) > IO_BKLASTEV(io)) {
# Get the next bucket.
if (qpio_rbucket (io, IO_EVI(io)) == EOF)
goto exit_
# Reject buckets that do not contain any events lying
# within the specified bounding box, if any.
if (bbused) {
ev_p = (IO_MINEVB(io) - 1) * SZ_SHORT / SZ_PIXEL + 1
$if (datatype == rd)
xs = Mem$t[ev_p+ev_xoff] + 0.5
ys = Mem$t[ev_p+ev_yoff] + 0.5
$else
xs = Mem$t[ev_p+ev_xoff]
ys = Mem$t[ev_p+ev_yoff]
$endif
ev_p = (IO_MAXEVB(io) - 1) * SZ_SHORT / SZ_PIXEL + 1
$if (datatype == rd)
xe = Mem$t[ev_p+ev_xoff] + 0.5
ye = Mem$t[ev_p+ev_yoff] + 0.5
$else
xe = Mem$t[ev_p+ev_xoff]
ye = Mem$t[ev_p+ev_yoff]
$endif
if (xs > x2 || xe < x1 || ys > y2 || ye < y1)
IO_EVI(io) = IO_BKLASTEV(io) + 1
}
}
# Copy out any events which pass the region mask and which share
# the same mask value. Note that in this case, to speed mask
# value lookup at random mask coordinates, the region mask for
# the bounding box is stored as a populated array in the QPIO
# descriptor.
ev = bp + (IO_EVI(io) - IO_BKFIRSTEV(io) - 1) * szs_event
bbmask = IO_BBMASK(io)
mval = 0
do i = IO_EVI(io), IO_BKLASTEV(io) {
# Get event x,y coordinates in whatever coord system.
ev = ev + szs_event
ev_p = (ev - 1) * SZ_SHORT / SZ_PIXEL + 1
$if (datatype == rd)
x = Mem$t[ev_p+ev_xoff] + 0.5
y = Mem$t[ev_p+ev_yoff] + 0.5
$else
x = Mem$t[ev_p+ev_xoff]
y = Mem$t[ev_p+ev_yoff]
$endif
# Reject events lying outside the bounding box.
if (bbused)
if (x < x1 || x > x2 || y < y1 || y > y2)
next
# Take a shortcut if no region mask is in effect for this BB.
if (bbmask == NULL)
goto putevent_
# Get the mask pixel associated with this event.
ii = (x - 1) / bb_xblock
jj = (y - 1) / bb_yblock
mval = Memi[bb_bufp + jj*bb_xsize + ii]
if (mval < 0)
mval = plr_getpix (bbmask, x, y)
# Accumulate points lying in the first nonzero mask range
# encountered.
if (mval != 0) {
if (maskval == 0)
maskval = mval
if (mval == maskval) {
putevent_ if (nev >= maxev)
break
nev = nev + 1
o_ev[nev] = ev
} else
break
}
}
IO_EVI(io) = i
case INDEX_NoRMorBB, INDEX_RMorBB:
# General extraction for indexed data. Process successive ranges
# and range lists until we get at least one event which lies within
# the bounding box, within a range, and which passes the event
# attribute filter, if one is in use.
# If the current range list (mask line) has been exhausted, advance
# to the next line which contains both ranges and events. A range
# list is used to specify the bounding box even if we don't have
# a nonempty region mask within the BB.
if (IO_RLI(io) > RLI_LEN(rl)) {
repeat {
y = IO_V(io,2)
if (IO_RLI(io) == RLI_INITIALIZE)
IO_RLI(io) = RL_FIRST
else
y = y + 1
if (y > IO_VE(io,2)) {
if (nev <= 0) {
o_nev = EOF
return (EOF)
} else
goto done_
}
IO_V(io,2) = y
evidx = Memi[IO_YOFFVP(io)+y-1]
if (evidx > 0) {
if (IO_RMUSED(io) == YES) {
if (IO_LINEIO(io) == YES) {
if (!pl_linenotempty (pl,IO_V(io,1)))
next
} else {
v[1] = IO_VE(io,1); v[2] = y
if (!pl_sectnotempty (pl,IO_V(io,1),v,NDIM))
next
}
call pl_glri (pl, IO_V(io,1), Memi[rl],
IO_MDEPTH(io), IO_VN(io,1), PIX_SRC)
}
IO_RLI(io) = RL_FIRST
}
} until (IO_RLI(io) <= RLI_LEN(rl))
IO_EVI(io) = evidx
IO_EV1(io) = evidx
IO_EV2(io) = Memi[IO_YLENVP(io)+y-1] + evidx - 1
}
# Refill the event bucket?
if (IO_EVI(io) > IO_BKLASTEV(io))
if (qpio_rbucket (io, IO_EVI(io)) == EOF)
goto exit_
# Compute current range parameters and initialize event pointer.
rp = rl + (IO_RLI(io) - 1) * RL_LENELEM
x1 = Memi[rp+RL_XOFF]
x2 = x1 + Memi[rp+RL_NOFF] - 1
maskval = Memi[rp+RL_VOFF]
ev = bp + (IO_EVI(io) - IO_BKFIRSTEV(io)) * szs_event
evtop = min (IO_EV2(io), IO_BKLASTEV(io))
# Extract events from bucket which lie within the current range
# of the current line. This is the inner loop of indexed event
# extraction, ignoring event attribute filtering.
do i = IO_EVI(io), evtop {
ev_p = (ev - 1) * SZ_SHORT / SZ_PIXEL + 1
$if (datatype == rd)
x = Mem$t[ev_p+ev_xoff] + 0.5
$else
x = Mem$t[ev_p+ev_xoff]
$endif
if (x >= x1) {
if (x > x2) {
IO_RLI(io) = IO_RLI(io) + 1
break
} else if (nev >= maxev)
break
nev = nev + 1
o_ev[nev] = ev
}
ev = ev + szs_event
}
IO_EVI(io) = i
if (i > IO_EV2(io))
IO_RLI(io) = RLI_NEXTLINE
}
done_
# Apply the event attribute filter if one is defined; repeat
# the whole process if we don't end up with any events.
if (nev > 0)
if (ex != NULL)
nev = qpex_evaluate (ex, o_ev, o_ev, nev)
if (nev <= 0)
goto again_
exit_
o_nev = nev
if (o_nev <= 0)
o_nev = EOF
return (o_nev)
end
$endfor
|