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
|
include "xyxymatch.h"
# RG_RDXYI -- Read in the x and y coordinates from a file and set the
# line number index.
int procedure rg_rdxyi (fd, x, y, lineno, xcolumn, ycolumn)
int fd #I the input file descriptor
pointer x #U pointer to the x coordinates
pointer y #U pointer to the y coordinates
pointer lineno #U pointer to the line numbers
int xcolumn #I column containing the x coordinate
int ycolumn #I column containing the y coordinate
int i, ip, bufsize, npts, lnpts, maxcols
pointer sp, str
real xval, yval
int fscan(), nscan(), ctor()
begin
call smark (sp)
call salloc (str, SZ_LINE, TY_CHAR)
bufsize = DEF_BUFSIZE
call malloc (x, bufsize, TY_REAL)
call malloc (y, bufsize, TY_REAL)
call malloc (lineno, bufsize, TY_INT)
maxcols = max (xcolumn, ycolumn)
npts = 0
lnpts = 0
while (fscan(fd) != EOF) {
lnpts = lnpts + 1
xval = INDEFR
yval = INDEFR
do i = 1, maxcols {
call gargwrd (Memc[str], SZ_LINE)
if (i != nscan())
break
ip = 1
if (i == xcolumn) {
if (ctor (Memc[str], ip, xval) <= 0)
xval = INDEFR
} else if (i == ycolumn) {
if (ctor (Memc[str], ip, yval) <= 0)
yval = INDEFR
}
}
if (IS_INDEFR(xval) || IS_INDEFR(yval))
next
Memr[x+npts] = xval
Memr[y+npts] = yval
Memi[lineno+npts] = lnpts
npts = npts + 1
if (npts >= bufsize) {
bufsize = bufsize + DEF_BUFSIZE
call realloc (x, bufsize, TY_REAL)
call realloc (y, bufsize, TY_REAL)
call realloc (lineno, bufsize, TY_INT)
}
}
call sfree (sp)
return (npts)
end
# RG_SORT -- If the coordinates are not already sorted, sort the coordinates
# first in y then in x. Remove points which are close together than a given
# tolerance, if the coincident point remove flag is on.
int procedure rg_sort (xcoord, ycoord, rsindex, npts, tolerance, sort, coincid)
real xcoord[ARB] #I pointer to the x coordinates
real ycoord[ARB] #I pointer to the y coordinates
int rsindex[ARB] #I pointer to sort index
int npts #I the number of objects
real tolerance #I coincidence tolerance in pixels
int sort #I sort the pixels ?
int coincid #I remove coincident points
int i, ndif
int rg_xycoincide()
begin
# Initialize the sort index.
do i = 1, npts
rsindex[i] = i
# Sort the pixels in y and then x if the arrays are unsorted.
if (sort == YES) {
call rg_qsortr (ycoord, rsindex, rsindex, npts)
call rg_sqsort (xcoord, ycoord, rsindex, rsindex, npts)
}
# Remove objects that are closer together than tolerance.
if (coincid == NO)
ndif = npts
else
ndif = rg_xycoincide (xcoord, ycoord, rsindex, rsindex, npts,
tolerance)
return (ndif)
end
# RG_XYCOINCIDE -- Remove points from a list which are closer together than
# a specified tolerance. The arrays are assumed to be sorted first in y then
# in x.
int procedure rg_xycoincide (xcoord, ycoord, a, b, npts, tolerance)
real xcoord[ARB] #I the input x coordinate values
real ycoord[ARB] #I the input y coordinate values
int a[ARB] #I the input sort index
int b[ARB] #O the output sort index
int npts #I the number of points
real tolerance #I the coincidence tolerace
int iprev, i, nunique
real tol2, r2
begin
tol2 = tolerance ** 2
nunique = npts
iprev = 1
repeat {
do i = iprev + 1, npts {
# Jump to the next object if this one has been deleted
# since all comparisons are then invalid.
if (a[iprev] == 0)
break
# Skip to the next object if this one has been deleted.
if (a[i] == 0)
next
# Check the tolerance limit in y and step to the next object
# if the bounds are exceeded.
r2 = (ycoord[a[i]] - ycoord[a[iprev]]) ** 2
if (r2 > tol2)
break
# Check the tolerance limit.
r2 = r2 + (xcoord[a[i]] - xcoord[a[iprev]]) ** 2
if (r2 <= tol2) {
a[i] = 0
nunique = nunique - 1
}
}
iprev = iprev + 1
} until (iprev >= npts)
# Reorder the index array.
if (nunique < npts) {
iprev = 0
do i = 1, npts {
if (a[i] != 0) {
iprev = iprev + 1
b[iprev] = a[i]
}
}
}
return (nunique)
end
|
