1 /* Copyright 2004 by Roger S. Bivand. */
2
3 #include "spdep.h"
4
polypoly(SEXP p1,SEXP n01,SEXP p2,SEXP n02,SEXP snap)5 SEXP polypoly(SEXP p1, SEXP n01, SEXP p2, SEXP n02, SEXP snap)
6 {
7 int n1=INTEGER_POINTER(n01)[0], n2=INTEGER_POINTER(n02)[0], pc=0;
8 int i, j, k=0;
9 double sn=NUMERIC_POINTER(snap)[0], dist;
10 double x1, x2, y1, y2, xd, yd;
11
12 SEXP ans;
13 PROTECT(ans = NEW_INTEGER(1)); pc++;
14
15 for (i=0; (i < n1) && (k < 2); i++) {
16 x1 = NUMERIC_POINTER(p1)[i];
17 y1 = NUMERIC_POINTER(p1)[n1 + i];
18 for (j=0; (j < n2) && (k < 2); j++) {
19 x2 = NUMERIC_POINTER(p2)[j];
20 y2 = NUMERIC_POINTER(p2)[n2 + j];
21 /* dist = pythag((x1-x2), (y1-y2));
22 if (dist < sn) k++;
23 if (k > 1) break;*/
24 /* following lines Micah Altman 2010 */
25 xd = x1-x2;
26 if (fabs(xd)>sn) { continue; }
27 yd = y1-y2;
28 if (fabs(yd)>sn) { continue; }
29 dist = hypot(xd, yd);
30 if (dist <= sn) {
31 k++;
32 }
33 }
34 }
35
36 INTEGER_POINTER(ans)[0] = k;
37
38 UNPROTECT(pc); /* ans */
39 return(ans);
40 }
41
42 /* function by Micah Altman */
43
spOverlap(SEXP bbbi,SEXP bbbj)44 SEXP spOverlap(SEXP bbbi, SEXP bbbj) {
45
46 int pc=0,overlap=1;
47 double bbi[4], bbj[4];
48 SEXP ans;
49
50 PROTECT(ans = NEW_INTEGER(1)); pc++;
51 bbi[0] = NUMERIC_POINTER(bbbi)[0];
52 bbi[1] = NUMERIC_POINTER(bbbi)[1];
53 bbi[2] = NUMERIC_POINTER(bbbi)[2];
54 bbi[3] = NUMERIC_POINTER(bbbi)[3];
55 bbj[0] = NUMERIC_POINTER(bbbj)[0];
56 bbj[1] = NUMERIC_POINTER(bbbj)[1];
57 bbj[2] = NUMERIC_POINTER(bbbj)[2];
58 bbj[3] = NUMERIC_POINTER(bbbj)[3];
59
60 if ((bbi[0]>bbj[2]) | (bbi[1]>bbj[3]) |
61 (bbi[2]<bbj[0]) | (bbi[3]<bbj[1]) ) {
62 overlap=0;
63 }
64
65 INTEGER_POINTER(ans)[0] = overlap;
66 UNPROTECT(pc); /* ans */
67 return(ans);
68 }
69
70 /* SEXP poly_loop(SEXP n, SEXP i_findInBox, SEXP bb, SEXP pl, SEXP nrs,
71 SEXP dsnap, SEXP criterion, SEXP scale) {
72
73 int nn = INTEGER_POINTER(n)[0];
74 int crit = INTEGER_POINTER(criterion)[0];
75 int Scale = INTEGER_POINTER(scale)[0];
76 int uBound = nn*Scale;
77 int i, j, jj, k, li, pc = 0;
78 int ii = 0;
79 int *card, *icard, *is, *jjs;
80
81 SEXP bbi, bbj, jhit, khit, ans, pli, plj, nrsi, nrsj;
82
83 int xx, yy, zz, ww;
84
85 card = (int *) R_alloc((size_t) nn, sizeof(int));
86 icard = (int *) R_alloc((size_t) nn, sizeof(int));
87 is = (int *) R_alloc((size_t) uBound, sizeof(int));
88 jjs = (int *) R_alloc((size_t) uBound, sizeof(int));
89
90 for (i=0; i<nn; i++) {
91 card[i] = 0;
92 icard[i] = 0;
93
94 }
95 for (i=0; i<uBound; i++) {
96 is[i] = 0;
97 jjs[i] = 0;
98 }
99
100 PROTECT(bbi = NEW_NUMERIC(4)); pc++;
101 PROTECT(bbj = NEW_NUMERIC(4)); pc++;
102 PROTECT(jhit = NEW_INTEGER(1)); pc++;
103 PROTECT(khit = NEW_INTEGER(1)); pc++;
104 PROTECT(nrsi = NEW_INTEGER(1)); pc++;
105 PROTECT(nrsj = NEW_INTEGER(1)); pc++;
106
107 for (i=0; i<(nn-1); i++) {
108 li = length(VECTOR_ELT(i_findInBox, i));
109 INTEGER_POINTER(nrsi)[0] = INTEGER_POINTER(nrs)[i];
110 for (k=0; k<4; k++)
111 NUMERIC_POINTER(bbi)[k] = NUMERIC_POINTER(bb)[i+(k*nn)];
112 for (j=0; j<li; j++) {
113 jj = INTEGER_POINTER(VECTOR_ELT(i_findInBox, i))[j] - ROFFSET;
114 for (k=0; k<4; k++)
115 NUMERIC_POINTER(bbj)[k] = NUMERIC_POINTER(bb)[jj+(k*nn)];
116 jhit = spOverlap(bbi, bbj);
117 if (INTEGER_POINTER(jhit)[0] > 0) {
118 INTEGER_POINTER(khit)[0] = 0;
119 INTEGER_POINTER(nrsj)[0] = INTEGER_POINTER(nrs)[jj];
120 if (INTEGER_POINTER(nrsi)[0]*INTEGER_POINTER(nrsj)[0] > 0){
121 khit = polypoly(VECTOR_ELT(pl, i), nrsi, VECTOR_ELT(pl, jj),
122 nrsj, dsnap);
123 }
124 if (INTEGER_POINTER(khit)[0] > crit) {
125 card[i]++;
126 card[jj]++;
127 is[ii] = i;
128 jjs[ii] = jj;
129 ii++;
130 if (ii == uBound) error("memory error, scale problem");
131 }
132 }
133 }
134 }
135
136 PROTECT(ans = NEW_LIST(nn)); pc++;
137
138 for (i=0; i<nn; i++) {
139 if (card[i] == 0) {
140 SET_VECTOR_ELT(ans, i, NEW_INTEGER(1));
141 INTEGER_POINTER(VECTOR_ELT(ans, i))[0] = 0;
142 } else {
143 SET_VECTOR_ELT(ans, i, NEW_INTEGER(card[i]));
144 }
145 }
146
147 for (i=0; i<ii; i++) {
148 xx = is[i];
149 yy = jjs[i];
150 zz = icard[yy];
151 ww = icard[xx];
152 if (zz == card[yy]) error("memory error, overflow");
153 if (ww == card[xx]) error("memory error, overflow");
154 INTEGER_POINTER(VECTOR_ELT(ans, yy))[zz] = xx + ROFFSET;
155 INTEGER_POINTER(VECTOR_ELT(ans, xx))[ww] = yy + ROFFSET;
156 icard[yy]++;
157 icard[xx]++;
158 }
159
160 for (i=0; i<nn; i++) {
161 if ((li = length(VECTOR_ELT(ans, i))) > 1) {
162 for (j=0; j<li; j++)
163 icard[j] = INTEGER_POINTER(VECTOR_ELT(ans, i))[j];
164 R_isort(icard, li);
165 for (j=0; j<li; j++)
166 INTEGER_POINTER(VECTOR_ELT(ans, i))[j] = icard[j];
167 }
168 }
169
170 UNPROTECT(pc);
171 return(ans);
172 } */
173
174
spOverlapC(double bbi1,double bbi2,double bbi3,double bbi4,double bbj1,double bbj2,double bbj3,double bbj4)175 int spOverlapC(double bbi1, double bbi2, double bbi3, double bbi4, double bbj1, double bbj2, double bbj3, double bbj4) {
176
177 int overlap=1;
178
179 if ((bbi1>bbj3) || (bbi2>bbj4) ||
180 (bbi3<bbj1) || (bbi4<bbj2) ) {
181 overlap = 0;
182 }
183
184 return(overlap);
185 }
186
polypolyC(double * px1,double * py1,int n1,double * px2,double * py2,int n2,double sn,int crit)187 int polypolyC(double *px1, double *py1, int n1, double *px2, double *py2,
188 int n2, double sn, int crit) {
189 int i, j, k=0;
190 double dist;
191 double x1, x2, y1, y2, xd, yd;
192
193 for (i=0; (i < n1) && (k < crit); i++) {
194 x1 = px1[i];
195 y1 = py1[i];
196 for (j=0; (j < n2) && (k < crit); j++) {
197 x2 = px2[j];
198 y2 = py2[j];
199 xd = x1-x2;
200 if (fabs(xd)>sn) { continue; }
201 yd = y1-y2;
202 if (fabs(yd)>sn) { continue; }
203 dist = hypot(xd, yd);
204 if (dist <= sn) {
205 k++;
206 }
207 }
208 }
209
210 return(k);
211 }
212
213
poly_loop2(SEXP n,SEXP i_findInBox,SEXP bb,SEXP pl,SEXP nrs,SEXP dsnap,SEXP criterion,SEXP nfIBB)214 SEXP poly_loop2(SEXP n, SEXP i_findInBox, SEXP bb, SEXP pl, SEXP nrs,
215 SEXP dsnap, SEXP criterion, SEXP nfIBB) {
216
217 int nn = INTEGER_POINTER(n)[0];
218 int crit = INTEGER_POINTER(criterion)[0];
219 /* int Scale = INTEGER_POINTER(scale)[0];*/
220 int uBound = (int) INTEGER_POINTER(nfIBB)[0]*2;
221 int i, j, jj, li, pc = 0;
222 int ii = 0;
223 int *card, *icard, *is, *jjs, *NRS, *cNRS;
224 double *bb1, *bb2, *bb3, *bb4, *plx, *ply;
225 double Dsnap = NUMERIC_POINTER(dsnap)[0];
226
227 // struct bbcontainer *bbs;
228
229 SEXP ans;
230
231 int jhit, khit, nrsi, nrsj;
232
233 int xx, yy, zz, ww;
234
235 card = (int *) R_alloc((size_t) nn, sizeof(int));
236 icard = (int *) R_alloc((size_t) nn, sizeof(int));
237 is = (int *) R_alloc((size_t) uBound, sizeof(int));
238 jjs = (int *) R_alloc((size_t) uBound, sizeof(int));
239 bb1 = (double *) R_alloc((size_t) nn, sizeof(double));
240 bb2 = (double *) R_alloc((size_t) nn, sizeof(double));
241 bb3 = (double *) R_alloc((size_t) nn, sizeof(double));
242 bb4 = (double *) R_alloc((size_t) nn, sizeof(double));
243 NRS = (int *) R_alloc((size_t) nn, sizeof(int));
244 cNRS = (int *) R_alloc((size_t) nn, sizeof(int));
245
246 for (i=0, li=0; i<nn; i++) {
247 card[i] = 0;
248 icard[i] = 0;
249 bb1[i] = NUMERIC_POINTER(bb)[i];
250 bb2[i] = NUMERIC_POINTER(bb)[i+(1*nn)];
251 bb3[i] = NUMERIC_POINTER(bb)[i+(2*nn)];
252 bb4[i] = NUMERIC_POINTER(bb)[i+(3*nn)];
253 NRS[i] = INTEGER_POINTER(nrs)[i];
254 li += NRS[i];
255 }
256
257 for (i=0; i<nn; i++) {
258 if (i == 0) cNRS[i] = 0;
259 else cNRS[i] = NRS[i-1] + cNRS[i-1];
260 }
261
262 for (i=0; i<uBound; i++) {
263 is[i] = 0;
264 jjs[i] = 0;
265 }
266
267 plx = (double *) R_alloc((size_t) li, sizeof(double));
268 ply = (double *) R_alloc((size_t) li, sizeof(double));
269
270 for (i=0, jj=0; i<nn; i++) {
271 nrsi = NRS[i];
272 for (j=0; j<nrsi; j++) {
273 plx[jj] = NUMERIC_POINTER(VECTOR_ELT(pl, i))[j];
274 ply[jj] = NUMERIC_POINTER(VECTOR_ELT(pl, i))[j+nrsi];
275 jj++;
276 /* if (i < (nn-1) && jj == li) error("polygon memory overflow");*/
277 }
278 }
279
280 for (i=0; i<(nn-1); i++) {
281 li = length(VECTOR_ELT(i_findInBox, i));
282 nrsi = NRS[i];
283 for (j=0; j<li; j++) {
284 jj = INTEGER_POINTER(VECTOR_ELT(i_findInBox, i))[j] - ROFFSET;
285 jhit = spOverlapC(bb1[i], bb2[i], bb3[i], bb4[i], bb1[jj],
286 bb2[jj], bb3[jj], bb4[jj]);
287 if (jhit > 0) {
288 khit = 0;
289 nrsj = NRS[jj];
290 if (nrsi > 0 && nrsj > 0){
291 khit = polypolyC(&plx[cNRS[i]], &ply[cNRS[i]], nrsi,
292 &plx[cNRS[jj]], &ply[cNRS[jj]], nrsj, Dsnap, crit+1L);
293 }
294 if (khit > crit) {
295 card[i]++;
296 card[jj]++;
297 is[ii] = i;
298 jjs[ii] = jj;
299 ii++;
300 /* if (ii == uBound) error("memory error, scale problem");*/
301 }
302 }
303 }
304 }
305
306 PROTECT(ans = NEW_LIST(nn)); pc++;
307
308 for (i=0; i<nn; i++) {
309 if (card[i] == 0) {
310 SET_VECTOR_ELT(ans, i, NEW_INTEGER(1));
311 INTEGER_POINTER(VECTOR_ELT(ans, i))[0] = 0;
312 } else {
313 SET_VECTOR_ELT(ans, i, NEW_INTEGER(card[i]));
314 }
315 }
316
317 for (i=0; i<ii; i++) {
318 xx = is[i];
319 yy = jjs[i];
320 zz = icard[yy];
321 ww = icard[xx];
322 /* if (zz == card[yy]) error("memory error, overflow");
323 if (ww == card[xx]) error("memory error, overflow");*/
324 INTEGER_POINTER(VECTOR_ELT(ans, yy))[zz] = xx + ROFFSET;
325 INTEGER_POINTER(VECTOR_ELT(ans, xx))[ww] = yy + ROFFSET;
326 icard[yy]++;
327 icard[xx]++;
328 }
329
330 for (i=0; i<nn; i++) {
331 if ((li = length(VECTOR_ELT(ans, i))) > 1) {
332 for (j=0; j<li; j++)
333 icard[j] = INTEGER_POINTER(VECTOR_ELT(ans, i))[j];
334 R_isort(icard, li);
335 for (j=0; j<li; j++)
336 INTEGER_POINTER(VECTOR_ELT(ans, i))[j] = icard[j];
337 }
338 }
339
340 UNPROTECT(pc);
341 return(ans);
342 }
343
344
345