1 /* Copyright where not otherwise stated Edzer Pebesma 2004,
2 copied from sp package */
3
4 # include <R.h>
5 # include <Rdefines.h>
6 /* # include <Rinternals.h> */
7 # define R_UNIFORM unif_rand()
8 # define R_NORMAL norm_rand()
9 # define RANDIN seed_in((long *) NULL)
10 # define RANDOUT seed_out((long *) NULL)
11
12 #ifndef MIN
13 # define MIN(a,b) ((a)>(b)?(b):(a))
14 #endif
15 #ifndef MAX
16 # define MAX(a,b) ((a)>(b)?(a):(b))
17 #endif
18
19 /* polygon structs: */
20 typedef struct {
21 double x, y;
22 } PLOT_POINT;
23
24 typedef struct {
25 PLOT_POINT min, max;
26 } MBR;
27
28 typedef struct polygon {
29 MBR mbr;
30 int lines;
31 PLOT_POINT *p;
32 int close; /* 1 - is closed polygon */
33 } POLYGON;
34
35 void setup_poly_minmax(POLYGON *pl);
36 static char InPoly(PLOT_POINT q, POLYGON *Poly);
37
R_point_in_polygon_mt(SEXP px,SEXP py,SEXP polx,SEXP poly)38 SEXP R_point_in_polygon_mt(SEXP px, SEXP py, SEXP polx, SEXP poly) {
39 int i/*, n*/;
40 PLOT_POINT p;
41 POLYGON pol;
42 SEXP ret;
43
44 pol.lines = LENGTH(polx); /* check later that first == last */
45 pol.p = (PLOT_POINT *) Calloc(pol.lines, PLOT_POINT); /* Calloc does error handling */
46 for (i = 0; i < LENGTH(polx); i++) {
47 pol.p[i].x = NUMERIC_POINTER(polx)[i];
48 pol.p[i].y = NUMERIC_POINTER(poly)[i];
49 }
50 pol.close = (pol.p[0].x == pol.p[pol.lines - 1].x &&
51 pol.p[0].y == pol.p[pol.lines - 1].y);
52 setup_poly_minmax(&pol);
53
54 ret = NEW_INTEGER(LENGTH(px));
55 for (i = 0; i < LENGTH(px); i++) {
56 p.x = NUMERIC_POINTER(px)[i];
57 p.y = NUMERIC_POINTER(py)[i];
58 /*
59 For each query point q, InPoly returns one of four char's:
60 i : q is strictly interior to P
61 o : q is strictly exterior to P
62 v : q is a vertex of P
63 e : q lies on the relative interior of an edge of P
64 */
65 switch (InPoly(p, &pol)) {
66 case 'i': INTEGER_POINTER(ret)[i] = 1; break;
67 case 'o': INTEGER_POINTER(ret)[i] = 0; break;
68 case 'v': INTEGER_POINTER(ret)[i] = 3; break;
69 case 'e': INTEGER_POINTER(ret)[i] = 2; break;
70 default: INTEGER_POINTER(ret)[i] = -1; break;
71 }
72 }
73 Free(pol.p);
74 return(ret);
75 }
76
setup_poly_minmax(POLYGON * pl)77 void setup_poly_minmax(POLYGON *pl) {
78 int i, n=pl->lines;
79 double minx,maxx,miny,maxy;
80
81 minx=miny=DBL_MAX;
82 maxx=maxy=-DBL_MAX;
83
84 for (i=0;i<n;i++) {
85 minx = MIN(minx, pl->p[i].x);
86 miny = MIN(miny, pl->p[i].y);
87 maxx = MAX(maxx, pl->p[i].x);
88 maxy = MAX(maxy, pl->p[i].y);
89 }
90 pl->mbr.min.x = minx;
91 pl->mbr.min.y = miny;
92 pl->mbr.max.x = maxx;
93 pl->mbr.max.y = maxy;
94 }
95
96 /*
97 This code is described in "Computational Geometry in C" (Second Edition),
98 Chapter 7. It is not written to be comprehensible without the
99 explanation in that book.
100
101 For each query point q, InPoly returns one of four char's:
102 i : q is strictly interior to P
103 o : q is strictly exterior to P
104 v : q is a vertex of P
105 e : q lies on the relative interior of an edge of P
106 These represent mutually exclusive categories.
107 For an explanation of the code, see Chapter 7 of
108 "Computational Geometry in C (Second Edition)."
109
110 Written by Joseph O'Rourke, contributions by Min Xu, June 1997.
111 Questions to orourke@cs.smith.edu.
112 --------------------------------------------------------------------
113 This code is Copyright 1998 by Joseph O'Rourke. It may be freely
114 redistributed in its entirety provided that this copyright notice is
115 not removed.
116 --------------------------------------------------------------------
117 */
118
119 /*
120 InPoly returns a char in {i,o,v,e}. See above for definitions.
121 */
122
InPoly(PLOT_POINT q,POLYGON * Poly)123 static char InPoly(PLOT_POINT q, POLYGON *Poly)
124 {
125 int n = Poly->lines;
126 PLOT_POINT *P=Poly->p;
127
128 int i, i1; /* point index; i1 = i-1 mod n */
129 double x; /* x intersection of e with ray */
130 double xx=q.x, yy=q.y;
131 int Rcross = 0; /* number of right edge/ray crossings */
132 int Lcross = 0; /* number of left edge/ray crossings */
133
134 /* For each edge e=(i-1,i), see if crosses ray. */
135 for( i = 0; i < n; i++ ) {
136 /* First see if q=(0,0) is a vertex. */
137 if (( P[i].x - xx )==0 &&( P[i].y - yy )==0 ) return 'v';
138 i1 = ( i + n - 1 ) % n;
139 /* printf("e=(%d,%d)\t", i1, i); */
140
141 /* if e "straddles" the x-axis... */
142 /* The commented-out statement is logically equivalent to the one
143 following. */
144 /* if( ( ( P[i].y > 0 ) && ( P[i1].y <= 0 ) ) ||
145 ( ( P[i1].y > 0 ) && ( P[i] .y <= 0 ) ) ) { }*/
146
147 if( (( P[i].y - yy ) > 0 ) != (( P[i1].y - yy ) > 0 ) ) {
148
149 /* e straddles ray, so compute intersection with ray. */
150 x = (( P[i].x - xx) *( P[i1].y - yy ) -( P[i1].x - xx ) *( P[i].y - yy )) /
151 (P[i1].y - P[i].y );
152 /* printf("straddles: x = %g\t", x); */
153
154 /* crosses ray if strictly positive intersection. */
155 if (x > 0) Rcross++;
156 }
157 /* printf("Right cross=%d\t", Rcross); */
158
159 /* if e straddles the x-axis when reversed... */
160 /* if( ( ( P[i] .y < 0 ) && ( P[i1].y >= 0 ) ) ||
161 ( ( P[i1].y < 0 ) && ( P[i] .y >= 0 ) ) ) { }*/
162
163 if ( (( P[i].y - yy ) < 0 ) != (( P[i1].y - yy ) < 0 ) ) {
164
165 /* e straddles ray, so compute intersection with ray. */
166 x = (( P[i].x - xx) *( P[i1].y - yy ) -( P[i1].x - xx ) *( P[i].y - yy )) /
167 (P[i1].y - P[i].y);
168 /* printf("straddles: x = %g\t", x); */
169
170 /* crosses ray if strictly positive intersection. */
171 if (x < 0) Lcross++;
172 }
173 /* printf("Left cross=%d\n", Lcross); */
174 }
175
176 /* q on the edge if left and right cross are not the same parity. */
177 if( ( Rcross % 2 ) != (Lcross % 2 ) )
178 return 'e';
179
180 /* q inside iff an odd number of crossings. */
181 if( (Rcross % 2) == 1 )
182 return 'i';
183 else return 'o';
184 }
185