1 | // Python - C extension for polygon module. |
---|
2 | // |
---|
3 | // To compile (Python2.3): |
---|
4 | // gcc -c polygon_ext.c -I/usr/include/python2.3 -o polygon_ext.o -Wall -O |
---|
5 | // gcc -shared polygon_ext.o -o polygon_ext.so |
---|
6 | // |
---|
7 | // See the module polygon.py |
---|
8 | // |
---|
9 | // |
---|
10 | // Ole Nielsen, GA 2004 |
---|
11 | // |
---|
12 | //NOTE: On 64 bit systems use long* instead of int* for Numeric arrays |
---|
13 | //this will also work on 32 bit systems |
---|
14 | |
---|
15 | |
---|
16 | #include "Python.h" |
---|
17 | #include "Numeric/arrayobject.h" |
---|
18 | #include "math.h" |
---|
19 | |
---|
20 | |
---|
21 | int _point_on_line(double x, double y, |
---|
22 | double x0, double y0, |
---|
23 | double x1, double y1) { |
---|
24 | /*Determine whether a point is on a line segment |
---|
25 | |
---|
26 | Input: x, y, x0, x0, x1, y1: where |
---|
27 | point is given by x, y |
---|
28 | line is given by (x0, y0) and (x1, y1) |
---|
29 | |
---|
30 | */ |
---|
31 | |
---|
32 | double a0, a1, a_normal0, a_normal1, b0, b1, len_a, len_b; |
---|
33 | |
---|
34 | a0 = x - x0; |
---|
35 | a1 = y - y0; |
---|
36 | |
---|
37 | a_normal0 = a1; |
---|
38 | a_normal1 = -a0; |
---|
39 | |
---|
40 | b0 = x1 - x0; |
---|
41 | b1 = y1 - y0; |
---|
42 | |
---|
43 | if ( a_normal0*b0 + a_normal1*b1 == 0 ) { |
---|
44 | //Point is somewhere on the infinite extension of the line |
---|
45 | |
---|
46 | len_a = sqrt(a0*a0 + a1*a1); |
---|
47 | len_b = sqrt(b0*b0 + b1*b1); |
---|
48 | |
---|
49 | if (a0*b0 + a1*b1 >= 0 && len_a <= len_b) { |
---|
50 | return 1; |
---|
51 | } else { |
---|
52 | return 0; |
---|
53 | } |
---|
54 | } else { |
---|
55 | return 0; |
---|
56 | } |
---|
57 | } |
---|
58 | |
---|
59 | |
---|
60 | |
---|
61 | int _separate_points_by_polygon(int M, // Number of points |
---|
62 | int N, // Number of polygon vertices |
---|
63 | double* points, |
---|
64 | double* polygon, |
---|
65 | long* indices, // M-Array for storage indices |
---|
66 | int closed, |
---|
67 | int verbose) { |
---|
68 | |
---|
69 | double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j; |
---|
70 | int i, j, k, outside_index, inside_index, inside; |
---|
71 | |
---|
72 | //Find min and max of poly used for optimisation when points |
---|
73 | //are far away from polygon |
---|
74 | |
---|
75 | minpx = polygon[0]; maxpx = minpx; |
---|
76 | minpy = polygon[1]; maxpy = minpy; |
---|
77 | |
---|
78 | for (i=0; i<N; i++) { |
---|
79 | px_i = polygon[2*i]; |
---|
80 | py_i = polygon[2*i + 1]; |
---|
81 | |
---|
82 | if (px_i < minpx) minpx = px_i; |
---|
83 | if (px_i > maxpx) maxpx = px_i; |
---|
84 | if (py_i < minpy) minpy = py_i; |
---|
85 | if (py_i > maxpy) maxpy = py_i; |
---|
86 | } |
---|
87 | |
---|
88 | //Begin main loop (for each point) |
---|
89 | inside_index = 0; //Keep track of points inside |
---|
90 | outside_index = M-1; //Keep track of points outside (starting from end) |
---|
91 | for (k=0; k<M; k++) { |
---|
92 | if (verbose){ |
---|
93 | if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M); |
---|
94 | } |
---|
95 | |
---|
96 | x = points[2*k]; |
---|
97 | y = points[2*k + 1]; |
---|
98 | |
---|
99 | inside = 0; |
---|
100 | |
---|
101 | //Optimisation |
---|
102 | if ((x > maxpx) || (x < minpx) || (y > maxpy) || (y < minpy)) { |
---|
103 | //Nothing |
---|
104 | } else { |
---|
105 | //Check polygon |
---|
106 | for (i=0; i<N; i++) { |
---|
107 | //printf("k,i=%d,%d\n", k, i); |
---|
108 | j = (i+1)%N; |
---|
109 | |
---|
110 | px_i = polygon[2*i]; |
---|
111 | py_i = polygon[2*i+1]; |
---|
112 | px_j = polygon[2*j]; |
---|
113 | py_j = polygon[2*j+1]; |
---|
114 | |
---|
115 | //Check for case where point is contained in line segment |
---|
116 | if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) { |
---|
117 | if (closed == 1) { |
---|
118 | inside = 1; |
---|
119 | } else { |
---|
120 | inside = 0; |
---|
121 | } |
---|
122 | break; |
---|
123 | } else { |
---|
124 | //Check if truly inside polygon |
---|
125 | if ( ((py_i < y) && (py_j >= y)) || |
---|
126 | ((py_j < y) && (py_i >= y)) ) { |
---|
127 | if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x) |
---|
128 | inside = 1-inside; |
---|
129 | } |
---|
130 | } |
---|
131 | } |
---|
132 | } |
---|
133 | if (inside == 1) { |
---|
134 | indices[inside_index] = k; |
---|
135 | inside_index += 1; |
---|
136 | } else { |
---|
137 | indices[outside_index] = k; |
---|
138 | outside_index -= 1; |
---|
139 | } |
---|
140 | } // End k |
---|
141 | |
---|
142 | return inside_index; |
---|
143 | } |
---|
144 | |
---|
145 | |
---|
146 | |
---|
147 | // Gateways to Python |
---|
148 | PyObject *point_on_line(PyObject *self, PyObject *args) { |
---|
149 | // |
---|
150 | // point_on_line(x, y, x0, y0, x1, y1) |
---|
151 | // |
---|
152 | |
---|
153 | double x, y, x0, y0, x1, y1; |
---|
154 | int res; |
---|
155 | PyObject *result; |
---|
156 | |
---|
157 | // Convert Python arguments to C |
---|
158 | if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1)) |
---|
159 | return NULL; |
---|
160 | |
---|
161 | |
---|
162 | // Call underlying routine |
---|
163 | res = _point_on_line(x, y, x0, y0, x1, y1); |
---|
164 | |
---|
165 | // Return values a and b |
---|
166 | result = Py_BuildValue("i", res); |
---|
167 | return result; |
---|
168 | } |
---|
169 | |
---|
170 | |
---|
171 | |
---|
172 | PyObject *separate_points_by_polygon(PyObject *self, PyObject *args) { |
---|
173 | //def separate_points_by_polygon(points, polygon, closed, verbose, one_point): |
---|
174 | // """Determine whether points are inside or outside a polygon |
---|
175 | // |
---|
176 | // Input: |
---|
177 | // point - Tuple of (x, y) coordinates, or list of tuples |
---|
178 | // polygon - list of vertices of polygon |
---|
179 | // closed - (optional) determine whether points on boundary should be |
---|
180 | // regarded as belonging to the polygon (closed = True) |
---|
181 | // or not (closed = False) |
---|
182 | |
---|
183 | // |
---|
184 | // Output: |
---|
185 | // indices: array of same length as points with indices of points falling |
---|
186 | // inside the polygon listed from the beginning and indices of points |
---|
187 | // falling outside listed from the end. |
---|
188 | // |
---|
189 | // count: count of points falling inside the polygon |
---|
190 | // |
---|
191 | // The indices of points inside are obtained as indices[:count] |
---|
192 | // The indices of points outside are obtained as indices[count:] |
---|
193 | // |
---|
194 | // Examples: |
---|
195 | // separate_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] ) |
---|
196 | // will return the indices [0, 2, 1] as only the first and the last point |
---|
197 | // is inside the unit square |
---|
198 | // |
---|
199 | // Remarks: |
---|
200 | // The vertices may be listed clockwise or counterclockwise and |
---|
201 | // the first point may optionally be repeated. |
---|
202 | // Polygons do not need to be convex. |
---|
203 | // Polygons can have holes in them and points inside a hole is |
---|
204 | // regarded as being outside the polygon. |
---|
205 | // |
---|
206 | // |
---|
207 | // Algorithm is based on work by Darel Finley, |
---|
208 | // http://www.alienryderflex.com/polygon/ |
---|
209 | // |
---|
210 | // |
---|
211 | |
---|
212 | PyArrayObject |
---|
213 | *points, |
---|
214 | *polygon, |
---|
215 | *indices; |
---|
216 | |
---|
217 | int closed, verbose; //Flags |
---|
218 | int count, M, N; |
---|
219 | |
---|
220 | // Convert Python arguments to C |
---|
221 | if (!PyArg_ParseTuple(args, "OOOii", |
---|
222 | &points, |
---|
223 | &polygon, |
---|
224 | &indices, |
---|
225 | &closed, |
---|
226 | &verbose)) |
---|
227 | return NULL; |
---|
228 | |
---|
229 | M = points -> dimensions[0]; //Number of points |
---|
230 | N = polygon -> dimensions[0]; //Number of vertices in polygon |
---|
231 | |
---|
232 | //FIXME (Ole): This isn't send to Python's sys.stdout |
---|
233 | if (verbose) printf("Got %d points and %d polygon vertices\n", M, N); |
---|
234 | |
---|
235 | //Call underlying routine |
---|
236 | count = _separate_points_by_polygon(M, N, |
---|
237 | (double*) points -> data, |
---|
238 | (double*) polygon -> data, |
---|
239 | (long*) indices -> data, |
---|
240 | closed, verbose); |
---|
241 | |
---|
242 | //NOTE: return number of points inside.. |
---|
243 | return Py_BuildValue("i", count); |
---|
244 | } |
---|
245 | |
---|
246 | |
---|
247 | |
---|
248 | // Method table for python module |
---|
249 | static struct PyMethodDef MethodTable[] = { |
---|
250 | /* The cast of the function is necessary since PyCFunction values |
---|
251 | * only take two PyObject* parameters, and rotate() takes |
---|
252 | * three. |
---|
253 | */ |
---|
254 | |
---|
255 | {"point_on_line", point_on_line, METH_VARARGS, "Print out"}, |
---|
256 | {"separate_points_by_polygon", separate_points_by_polygon, |
---|
257 | METH_VARARGS, "Print out"}, |
---|
258 | {NULL, NULL, 0, NULL} /* sentinel */ |
---|
259 | }; |
---|
260 | |
---|
261 | |
---|
262 | |
---|
263 | // Module initialisation |
---|
264 | void initpolygon_ext(void){ |
---|
265 | Py_InitModule("polygon_ext", MethodTable); |
---|
266 | |
---|
267 | import_array(); //Necessary for handling of NumPY structures |
---|
268 | } |
---|
269 | |
---|
270 | |
---|
271 | |
---|