1 # include <R.h>
2 # include <Rdefines.h>
3
4 #define ROFFSET 1
5
6 int pipbb(double pt1, double pt2, double *bbs);
7
8 int between(double x, double low, double up);
9
mtInsiders(SEXP n1,SEXP bbs)10 SEXP mtInsiders(SEXP n1, SEXP bbs) {
11
12 int n, pc=0;
13 int i, j, k, k1;
14 double bbi[4], bbj[4];
15 int *yes, jhit[4], hsum;
16 SEXP ans;
17
18 n = INTEGER_POINTER(n1)[0];
19 PROTECT(ans = NEW_LIST(n)); pc++;
20 yes = (int *) R_alloc((size_t) n, sizeof(int));
21 for (i=0; i < n; i++) {
22 for (j=0; j < n; j++) yes[j] = 0;
23 bbi[0] = NUMERIC_POINTER(bbs)[i];
24 bbi[1] = NUMERIC_POINTER(bbs)[i + n];
25 bbi[2] = NUMERIC_POINTER(bbs)[i + 2*n];
26 bbi[3] = NUMERIC_POINTER(bbs)[i + 3*n];
27 k = 0;
28 for (j=0; j < n; j++) {
29 if (i != j) {
30 hsum = 0;
31 bbj[0] = NUMERIC_POINTER(bbs)[j];
32 bbj[1] = NUMERIC_POINTER(bbs)[j + n];
33 bbj[2] = NUMERIC_POINTER(bbs)[j + 2*n];
34 bbj[3] = NUMERIC_POINTER(bbs)[j + 3*n];
35 for (k1=0; k1 < 4; k1++) jhit[k1] = 0;
36 jhit[0] = pipbb(bbi[2], bbi[3], bbj);
37 jhit[1] = pipbb(bbi[0], bbi[1], bbj);
38 jhit[2] = pipbb(bbi[0], bbi[3], bbj);
39 jhit[3] = pipbb(bbi[2], bbi[1], bbj);
40
41 for (k1=0; k1 < 4; k1++)
42 hsum = hsum + jhit[k1];
43 if (hsum == 4) { yes[j] = 1;
44 k = k + yes[j];
45 }
46 }
47 }
48
49 if (k != 0) {
50 SET_VECTOR_ELT(ans, i, NEW_INTEGER(k));
51 for (j=0, k1=0; j < n; j++) {
52 if (yes[j] > 0)
53 INTEGER_POINTER(VECTOR_ELT(ans,
54 i))[k1++] = j + ROFFSET;
55 }
56 }
57 }
58 UNPROTECT(pc); /* ans */
59 return(ans);
60 }
61
62
between(double x,double low,double up)63 int between(double x, double low, double up) {
64 if (x >= low && x <= up) return(1);
65 else return(0);
66 }
67
pipbb(double pt1,double pt2,double * bbs)68 int pipbb(double pt1, double pt2, double *bbs) {
69 if ((between(pt1, bbs[0], bbs[2]) == 1) &&
70 (between(pt2, bbs[1], bbs[3]) == 1)) return(1);
71 else return(0);
72 }
73
74